Phase transitions in self-gravitating systems and bacterial populations 
with a screened attractive potential 
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We consider a system of particles interacting via a screened Newtonian potential and study phase 
transitions between homogeneous and inhomogeneous states in the microcanonical and canonical 
ensembles. Like for other systems with long-range interactions, we obtain a great diversity of 
microcanonical and canonical phase transitions depending on the dimension of space and on the 
importance of the screening length. We also consider a system of particles in Newtonian interaction 
in the presence of a "neutralizing background" . By a proper interpretation of the parameters, our 
study describes (i) self-gravitating systems in a cosmological setting, and (ii) chemotaxis of bacterial 
populations in the original Keller-Segel model. 
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I. INTRODUCTION 



Many biological species like bacteria, amoebae, en- 
dothelial cells, or even ants interact through the phe- 
nomenon of chemotaxis pj. These organisms secrete a 
chemical substance (like a pheromone) that has an at- 
tractive (or sometimes repulsive) action on the organ- 
isms themselves. This phenomenon is responsible for the 
self-organization and morphogenesis of many biological 
species. It has also been proposed as a leading mechanism 
for the formation of blood vessels during embriogenesis 
• On a theoretical point of view, chemotaxis can be de- 
scribed by the Keller-Segel [|| model or its generalizations 
[H. The Keller-Segel model consists in a drift-diffusion 
equation for the evolution of the density of bacteria p(r, t) 
coupled to a reaction-diffusion equation for the evolu- 
tion of the secreted chemical c(r,t). In certain approxi- 
mations, the reaction-diffusion equation is replaced by a 
Poisson equation. In that case, the Keller-Segel (KS) 
model becomes isomorphic to the Smoluchowski-Poisson 
(SP) system Q describing self-gravitating Brownian par- 
ticles (see, e.g., [f| for a description of this analogy). The 
KS model and SP system have been studied thoroughly 
in applied mathematics (see Rcfs. in 0) and in theoret- 
ical physics (see Refs. in [Bj]). 

However, the original KS model Q also allows for the 
possibility that the chemical suffers a degradation process 
which has the effect of reducing the range of the interac- 
tion. In that case, the Poisson equation is replaced by a 
screened Poisson equation Q. In the gravitational anal- 
ogy, this amounts to replacing the gravitational potential 
by a screened gravitational potential, i.e. an attractive 
Yukawa potential. In that case, there exists interest- 
ing phase transitions between spatially homogeneous and 
spatially inhomogeneous equilibrium distributions. This 
is a physical motivation to consider the thermodynamics 
of N-body systems interacting via an attractive Yukawa 
potential (9f. This will be called the screened Newto- 
nian model. We shall also consider a related model where 
the interaction is not screened but the Poisson equation 



is modified so as to allow for the existence of spatially 
homogeneous distributions at equilibrium. This will be 
called the modified Newtonian model. In that model, the 
source of the potential is the deviation between the ac- 
tual density p(r, t) and the average density p. This is 
similar to the effect of a "neutralizing background" in 
plasma physics . This model can be derived from the 
Keller-Segel model in the limit of vanishing degradation 
of the chemical [ll|. It also appears in cosmology, due 
to the expansion of the universe, when we work in the 
comoving frame [Hj]. It is therefore interesting to con- 
sider this form of interaction at a general level and study 
the corresponding phase transitions. We shall also com- 
pare them with the ones obtained within the ordinary 
Newtonian model [l3l - l30j (see a review in (3l| ). 

The paper is organized as follows. In Sec. UH we 
discuss several kinetic models taken from astrophysics, 
plasma physics and biology for which our study applies. 
We consider either isolated systems described by the mi- 
crocanonical ensemble (fixed energy E) or dissipative sys- 
tems described by the canonical ensemble (fixed temper- 
ature T). We characterize their equilibrium states in the 
mean field approximation: in the microcanonical ensem- 
ble (MCE), they maximize the entropy at fixed mass and 
energy and in the canonical ensemble (CE) they min- 
imize the free energy at fixed mass. In Sec. IIII1 we 
specifically consider the case of a Newtonian interaction 
with a neutralizing background. We study phase tran- 
sitions between homogeneous and inhomogeneous states 
depending on the dimension of space. In d = 1, the 
system presents canonical and microcanonical second or- 
der phase transitions. In d = 2, the system presents 
an isothermal collapse in CE (zeroth order phase transi- 
tion) and a first order phase transition in MCE. In d = 3, 
the system presents an isothermal collapse in CE and a 
gravothermal catastrophe in MCE (zeroth order phase 
transitions). In Sec. IIV1 we perform a similar study 
for the attractive Yukawa potential with screening length 
Icq . In d = 1, there exists a canonical tricritical point 
(ko) c R = V^tt ~ 4.44 and a microcanonical tricritical 
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point (ko) m R ~ 11.8, where R is the system size. If 
< (fco)cj the system presents canonical and micro- 
canonical second order phase transitions. In that case, 
the ensembles are equivalent. If (fco) c < h) < (ko) m , 
the system presents a canonical first order phase transi- 
tion and a microcanonical second order phase transition. 
In that case, there exists a region of negative specific 
heats in MCE and the ensembles are inequivalcnt. If 
kg > (fco)roi the system presents canonical and micro- 
canonical first order phase transitions. In d = 2 and 
d = 3, the phase transitions are similar to those reported 
for the modified Newtonian model. In Sec. [V] we study 
the dynamical stability of the homogeneous phase and 
analytically determine the critical point (E*,T*) that 
marks the onset of instability of the homogeneous branch 
and the starting point of the bifurcated inhomogeneous 
branch. Direct numerical simulations associated with 
these phase transitions will be reported in a forthcom- 
ing paper. 

Finally, it may be noted that the phase transitions 
reported in this paper share analogies (but also differ- 
ences) with phase transitions observed in the Hamil- 
tonian mean field (HMF) model f32^36] . the spherical 
mass shell (SMS) model (37|, the Blume-Emery-Griffiths 
(BEG) model 13811 . the infinite-range attactive interaction 
(IRAI) model [39||, the self-gravitating Fermi gas (SGF) 
model the self-gravitating ring (SGR) model |40j 
and the one-dimensional static cosmology (OSC) model 

m. 



with 



S = —kn I — In — drdv 



M = I pdr, 



(4) 



(5) 



E = J /^drrfv+i J p(v,t)u(v,v')p(r',t) dvdv' 

pVdr, (6) 



where p(r,t) = J f(r,v,t)dv is the spatial density. In- 
troducing the mean field potential 

$(r) = J u(r,r')p(r')dr' + V(r), (7) 

the energy can also be written 

E = - J fv 2 drdv + i J p($ + V) dr. (8) 

We shall be interested in global and local entropy max- 
ima. Let us first determine the critical points of entropy 
at fixed mass and energy which cancel the first order vari- 
ations. Introducing Lagrange multipliers, they satisfy 



II. KINETIC MODELS AND STATISTICAL 
EQUILIBRIUM STATES 

1, Isolated systems 



SS - -5E - aSM = 0. 



(9) 



The variations are straightforward to evaluate and we 
obtain the mean field Maxwcll-Boltzmann distribution 



We consider an isolated system of N particles in inter- 
action described by the Hamiltonian equations 



dri dH dvi dH 
dt dvj ' dt dri ' 



(1) 



where 

i i<j i 

We assume that the particles interact through a binary 
potential u(r, r') that is symmetric with respect to the 
interchange of r and r', and that they also evolve in a 
fixed external potential V(r). Since the system is iso- 
lated, with strict conservation of energy and mass, the 
proper statistical ensemble is the microcanonical ensem- 
ble Q. In this paper, we shall use a mean field approach 
[63| . In the microcanonical ensemble, the statistical equi- 
librium state is obtained by maximizing the Boltzmann 
entropy at fixed mass and energy. We thus have to solve 
the maximization problem 



max{5[/] 



E[f]=E,M[f]=M}, 



(3) 



/ = Ae 



(10) 



where j3 = 1/ksT and $(r) is given by Eq. J7]). Integrat- 
ing over the velocity, the find that the density is given by 
the mean field Boltzmann distribution 



p = A'e 



(11) 



This critical point is a (local) entropy maximum at fixed 
mass and energy iff 



S 2 J = - 



^X- drd v - -P I 8p5<P dr < 0, (12) 
2mf 2 



for all perturbations 5f that conserve mass and energy 
at first order. In Appendix \K\ we provide an equivalent 
but simpler condition of stability in the microcanonical 
ensemble [see inequality (|A12I) ]. 

The time evolution of the distribution function 
/(r, v, t) is governed by a kinetic equation of the form 



d -L + v d -i-v*. d -l = ( d -l 

dt dv dr V dt 



(13) 
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where 



with 



$(r,t) - J u(r,r')p(r',t) dr' + V(r), 



(14) 



is the time-dependent mean field potential. The l.h.s. is 
an advective operator (Vlasov) in phase space. The r.h.s. 
is a "collision" operator like the Boltzmann operator in 
the kinetic theory of gases or like the Landau (or Lenard- 
Balcscu) operator in plasma physics or stellar dynamics. 
The "collision" operator in Eq. (fT3")) takes into account 
the development of correlations between particles. It 
can have a more or less complicated form but it satis- 
fies general properties associated with the first and sec- 
ond principles of thermodynamics: (i) it conserves mass 
and energy; (ii) it satisfies an iJ-theorem for the Boltz- 
mann entropy (j4|), i.e. S > with an equality iff / is 
the Maxwell-Boltzmann distribution (|10p . Furthermore, 
the Maxwell-Boltzmann distribution is dynamically sta- 
ble iff it is a (local) entropy maximum at fixed mass and 
energy. These general properties can be checked directly 
for the Boltzmann equation, for the Landau equation, for 
the Lenard-Balescu equation and for the BGK operator. 
Therefore, the kinetic equation (|13[) is consistent with 
the maximization problem ([3]) describing the statistical 
equilibrium state of the system in MCE. If we neglect the 
collisions for sufficiently short times, Eq. (fT3"f reduces to 
the Vlasov equation which can experience a complicated 
process of collisionless violent relaxation towards a quasi 
stationary state (QSS) [43j. 



2. Dissipative systems in phase space 

We consider a dissipative system of N Brownian par- 
ticles in interaction described by the Langevin equations 



dvj 
dt 



dri 

m 

dt 

1 dH 
m dr i 



dH 



^Ri(t), 



(15) 



(16) 



where H is the Hamiltonian defined by Eq. ([2]), — £v, 
is a friction force and Rj(i) is a white noise satisfying 
(Ri(t)) = and {R^{t)R^(t)} = S tj 6^5(t-t'). The diffu- 
sion coefficient D and the friction coefficient £ are related 
to each other according to the Einstein relation £ = Dj3m 
where /3 = l/(fcsT) is the inverse temperature. Since 
this system is dissipative, the proper statistical ensemble 
is the canonical ensemble 0] . In the canonical ensemble, 
the statistical equilibrium state is obtained by minimiz- 
ing the Boltzmann free energy F[f] = E[f] — TS[f] at 
fixed mass. We thus have to solve the minimization prob- 
lem 



nun 



M[f] = M} 



(17) 



F= //ydrdv+I J p(r,t)u(v,v')p(v',t)drdr' 

+ f pV dr + k B T f — In — drdv. (18) 
J J m m 

We shall be interested by global and local minima of free 
energy. Let us first determine the critical points of free 
energy at fixed mass which cancel the first order varia- 
tions. Introducing a Lagrange multiplier, they satisfy 



5F + aTSM = 0. 



(19) 



The variations are straightforward to evaluate and we 
obtain the mean field Maxwell-Boltzmann distribution 
(p~0|l and the mean field Boltzmann distribution (TlTT) as 
in the microcanonical ensemble. This critical point is a 
(local) minimum of free energy iff 



S 2 F 



1 



Sp8& dr 



m 

2/ 



drdv > 0, (20) 



for all perturbations df that conserve mass. 

In the mean field approximation, the evolution of the 
distribution function f(r,v,t) is governed by a kinetic 
equation of the form 



df df „^ df d f df 

— + v • — - V$ ■ — = — D— + £/v 

dt dv dr dv V dv 



(21) 



coupled to the mean field potential (fbfj). This is called 
the mean field Kramers equation. The mean field 
Kramers equation conserves mass and satisfies an H- 
theorem for the Boltzmann free energy (TT5)) . i.e. F < 
with an equality iff / is the Maxwell-Boltzmann distribu- 
tion (pTO]) . Furthermore, the Maxwell-Boltzmann distri- 
bution is dynamically stable iff it is a (local) minimum of 
free energy at fixed mass. Therefore, the kinetic equation 
(|2ip is consistent with the minimization problem (JTTJ) de- 
scribing the statistical equilibrium state of the system in 
CE. 

Remark: the critical points in MCE and CE are the 
same because the variational problems ^ and ([TT]) arc 
equivalent at the level of the first order variations ([9]) and 
(Unj). However, they are not equivalent at the level of the 
second order variations (fT2"j) and (|2H|) because of the dif- 
ferent class of perturbations to consider. Therefore, we 
can have ensembles inequivalence [22l. l3ll. [44l |45| . In fact, 
the condition of canonical stability (|17[) provides a suf- 
ficient condition of microcanonical stability @. Indeed, 
if inequality (|20l) is satisfied for all perturbations that 
conserve mass, then it is a fortiori satisfied for perturba- 
tions that conserve mass and energy, so that inequality 
(fT2"| is satisfied. Therefore, canonical stability implies 
microcanonical stability: 



CE3) => ©. 



(22) 



However, the converse is wrong in case of ensembles in- 
equivalence. 
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3. Dissipative systems in physical space 

In the strong friction limit £ — > +00, we can formally 
neglect the inertial term dvi /dt in Eq. (fTB"]) and we obtain 
the overdampcd Langcvin equations 



~dt 



/ 2DTL,{t). 



(23) 



The statistical equilibrium state of this system (described 
by the canonical ensemble ) is obtained by solving the 
minimization problem 



with 



mm{F[p]\M[p] = M} , 
p 



F =\ j p(r,t)u(r,r')p(r',t)drdr' 
pVdr + k B T [ — In — dr. 



(24) 



Writing the variational principle as 
SF + aTSM = 0, 



(25) 



(26) 



we obtain the mean field Boltzmann distribution (1111) . 
This critical point is a (local) minimum of free energy at 
fixed mass iff 



52p = \ I 



k B T f(Sp) 



2p 



dr > 0, (27) 



for all perturbations dp that conserves mass. 

In the mean field approximation, the evolution of the 
density profile p(r, t) is governed by a kinetic equation of 
the form 



dp 
dt 



k R T. 



-V/9 + pV$ 



(28) 



coupled to the mean field equation (fT4"|) . This is called 
the mean field Smoluchowski equation. The mean field 
Smoluchowski equation (|28|) conserves mass and satisfies 
an _H-thcorcm for the Boltzmann free energy (|25[) . i.e. 
F < with an equality iff p is the Boltzmann distribu- 
tion (|11|). Furthermore, the Boltzmann distribution is 
dynamically stable iff it is a (local) minimum of free en- 
ergy at fixed mass. Therefore, the kinetic equation ([28]) is 
consistent with the minimization problem (|24[) describing 
the statistical equilibrium state of the system in CE. 

Remark 1: the Smoluchowski equation (|28[) can also be 
deduced from the Kramers equation (|21[) in the strong 
friction limit [46|. For £, D —t +00 and /3 = £/Dm fi- 
nite, the time-dependent distribution function /(r, v, t) 
is Maxwellian 



/(r,v,t)= ( — 



d/2 



p(r,t)e 



0(1/0, (29) 



and the time-dependent density p(r, t) is solution of the 
Smoluchowski equation (f2"5]l . Using Eq. ([2U]) , we can ex- 
press the free energy (fT5| as a functional of the density 
and we obtain the free energy (|25p up to some unimpor- 
tant constants. 

Remark 2: it is shown in Appendix [X] that the max- 
imization problems (|17|) and (|24[) arc equivalent in the 
sense that /(r, v) is solution of (fT?]) iff p(r) is solution of 
(1241). Thus, we have 



02) <^ il. 



(30) 



As a consequence, the Maxwcll-Boltzmann distribution 
/(r,v) is dynamically stable with respect to the mean 
field Kramers equation ([2Tj) iff the corresponding Boltz- 
mann distribution p(r) is dynamically stable with re- 
spect to the mean field Smoluchowski equation (|28|) . On 
the other hand, according to the implication (f22"|) . the 
Maxwcll-Boltzmann distribution /(r,v) is dynamically 
stable with respect to the kinetic equation (TT3]) if it is 
stable with respect to the mean field Kramers equation 
([2"T]) . but the reciprocal is wrong in case of ensembles 
inequivalence. 



4- The Keller-Segel model of chemotaxis 

The Keller-Segel model [|| describing the chemotaxis 
of biological populations can be written as 



dp 

dt 



V ■ (DVp - xpVc) 



\_dc 
LTdt 



Ac - k 2 c + Xp, 



(31) 



(32) 



where p is the concentration of the biological species (e.g. 
bacteria) and c is the concentration of the secreted chem- 
ical. The bacteria diffuse with a diffusion coefficient D 
and undergo a chemotactic drift with strength \ along 
the gradient of chemical. The chemical is produced by 
the bacteria at a rate D'X, is degraded at a rate D'k 2 
and diffuses with a diffusion coefficient D'. We adopt 
Neumann boundary conditions Q: 



Vc ■ n = 0, Vp ■ n = 0, 



(33) 



where n is a unit vector normal to the boundary of the 
domain. The drift-diffusion equation pip is similar to 
the mean field Smoluchowski equation ([28]) where the 
concentration of chemical — c(r, t) plays the role of the 
potential $(r, t). Therefore, there exists many analogies 
between chemotaxis and Brownian particles in interac- 
tion 0. In particular, the effective statistical ensemble 
associated with the Keller-Segel model is the canonical 
ensemble. The steady states of the Keller-Segel model 
are of the form 



p = Ae~n l 



(34) 
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which is similar to the Boltzmann distribution (fTTj) with 
an effective temperature T e ff = D/x- The Lyapunov 
functional associated with the KS model is Q : 

F ^i\I f (Vc)2 + k2c ^ dr ~ I pcdr 

+T eff J plnpdr. (35) 

It is similar to a free energy F = E — T e ffS in thermo- 
dynamics, where E is the energy and S is the Boltzmann 
entropy The KS model conserves mass and satisfies an 
-ff-theorem for the free energy (|35p . i.e. F < with an 
equality iff p is the Boltzmann distribution (f3~4"|) . Further- 
more, the Boltzmann distribution is dynamically stable 
iff it is a (local) minimum of free energy at fixed mass. 
In that context, the minimization problem 



min {F[p, c] 

p,c 



M[p] = M\ . 



(36) 



determines a steady state of the KS model that is dy- 
namically stable. This is similar to a condition of ther- 
modynamical stability in the canonical ensemble. 

Let us consider some simplified forms of the Keller- 
Segel model that have been introduced in the literature: 
(i) In the limit of large diffusivity of the chemical D' —> 
+oo at fixed k and A, the reaction-diffusion equation 
takes the form of a screened Poisson equation Q: 



Ac 



k 2 c 



-Xp, 



(37) 



and the free energy becomes 

F = --J pcdr + T eff J plnpdr. (38) 

In that case, the KS model is isomorphic to the Smolu- 
chowski equation (|28[) with an attractive Yukawa poten- 
tial flu}- 

(ii) In the limit of large diffusivity of the chemical 
D' — > +oo and a vanishing degradation rate k 2 = 0, 
the reaction-diffusion equation (|32l) takes the form of a 
modified Poisson equation : 



Ac=-X(p-p), 



(39) 



where p = M/V is the average density, and the free en- 
ergy becomes 

F =-\j(p- P)cdr + T eff J plnpdr. (40) 

In that case, the KS model is isomorphic to the Smolu- 
chowski equation (|28[) with a modified Poisson equation 
©). 

(iii) Some authors have also considered a simple model 
of chemotaxis where the reaction-diffusion equation (|32[) 
is replaced by the Poisson equation [47J : 



This is valid in the absence of degradation of the chemical 
and for sufficiently large densities p> p. This model can 
be used in particular to study chemotactic collapse. The 
corresponding free energy is 



F 



- I pcdr + T eff I plnpdr. 



(42) 



Ac 



-Xp. 



(41) 



In that model, the boundary conditions (|33[) must be 
modified [iH and we must impose that c — ► at infinity 
like for the gravitational potential in astrophysics. Fur- 
thermore, we must impose that the normal component of 
the current vanishes on the boundary: (Z?V p — XP^ C ) ' 
n = so as to conserve mass. In that case, the KS 
model is isomorphic to the Smoluchowski-Poisson (SP) 
system describing self-gravitating Brownian particles in 
the overdamped limit [29(. 



5. Physical justification of the canonical ensemble for 
systems with long-range interactions 

In statistical mechanics, the canonical distribution is 
usually derived by considering a subpart of a large sys- 
tem and assuming that the rest of the system plays the 
role of a thermostat (H)]. However, this justification im- 
plicitly assumes that energy is additive. Since energy is 
non-additive for systems with long-range interactions, it 
is sometimes concluded that the canonical ensemble has 
no foundation to describe systems with long-range inter- 
actions [49| . In fact, this is not quite true @. We can give 
two justifications of the canonical ensemble for systems 
with long-range interactions: 

(i) The canonical ensemble is relevant to describe a 
system of particles in contact with a thermal bath of a 
different nature 0. This is the case if we consider a 
system of Brownian particles in interaction described by 
the stochastic equations (|T5|l - (ri7J)l . The particles interact 
through a potential u(r, r') that can be long-range, but 
they also undergo a friction force and a stochastic force 
that are due to other types of interaction (they model 
in general short-range interactions). As we have seen, 
this system is described by the canonical ensemble. It 
does not correspond to a subsystem of a larger system, 
but simply to a system as a whole with long-range and 
short-range interactions (65[. 

(ii) Since canonical stability implies microcanonical 
stability [HJ , the condition of canonical stability provides 
a sufficient condition of microcanonical stability. In this 
sense, the canonical stability criterion (see Sees. Ill 2 1 and 
III 3[) can be useful even for an isolated Hamiltonian sys- 
tem (see Sec. Ill 1|) because if we can prove that this 
system is canonically stable, then it is granted to be mi- 
crocanonically stable. This remark also applies to other 
ensembles (grand canonical, grand microcanonical,...). 
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III. THE MODIFIED NEWTONIAN MODEL 

In this section, we discuss phase transitions that ap- 
pear in the modified Newtonian model. 

A. Physical motivation of the model 

We consider a system of particles interacting via a 
mean field potential $(r, f) that is solution of the modi- 
fied Poisson equation 



A$ = S d G( P -p), 



(43) 



where p = M/V is the average density (conserved quan- 
tity). At statistical equilibrium, the density is given by 
the Boltzmann distribution 



(44) 



We have used the notations of astrophysics (where G is 
the constant of gravity and the surface of a unit sphere 
in d dimensions) in order to make the connection with 
ordinary self-gravitating systems where Eq. (|43[) is re- 
placed by the Poisson equation A<£> = SdGp. However, 
this model can have application in other contexts as ex- 
plained below. We assume that the system is confined 
in a finite domain (box) and we impose the Neumann 
boundary conditions 



V$-n = 0, Vp-n = 0, 



(45) 



where n is a unit vector normal to the boundary of the 
box (the explicit expression of the potential in d = 1 is 
given in Appendix |B|) . This model admits spatially ho- 
mogeneous solutions (p = ~p and $ = 0) at any tempera- 
ture. It also admits spatially inhomogeneous solutions at 
sufficiently low temperatures. We shall study this model 
in arbitrary dimensions of space d with explicit compu- 
tations for d = 1, 2, 3. This model has different physical 
applications: 

(i) It describes self-gravitating systems in a cosmolog- 
ical setting 12j. Due to the expansion of the universe, 
when we work in the comoving frame, the Poisson equa- 
tion takes the form of Eq. (|43|) where the potential is 
produced by the deviation between the actual density 
p(r,t) and the mean density p. In cosmology, we must 
also account for the scale factor a(t) but if we consider 
timcscalcs that are short with respect to the Hubble time 
H^ 1 = a/ a, we can ignore this time dependence. This 
model has been studied by Valageas [411 in d — 1 with pe- 
riodic boundary conditions. In that context, the relevant 
ensemble is the MCE since the system is isolated. 

(ii) By a proper reinterpretation of the parameters, 
the field equation (|4*3"]) describes the relation between the 
concentration of the chemical and the density of bacteria 
in the Kcllcr-Segel model (f3T)l) . In that case, the most 
physical dimension is d — 2 and the boundary conditions 
are of the form (|45[) . Furthermore, the relevant ensemble 



is the CE since the KS model has a canonical structure. 
This model has been studied by applied mathematicians, 
starting with Jagcr & Luckhaus PJJ, but they have not 
performed the type of study that we are developing in 
this paper. 

In view of these different applications, we shall study 
this model in the microcanonical and canonical ensembles 
in any dimension of space. 



B. The modified Emden equation 

In the modified Newtonian model, the statistical equi- 
librium state is given by the Boltzmann distribution (|44[) 
coupled to the modified Poisson equation (|4*3")l . We look 
for spherically symmetric solutions because, for non ro- 
tating systems, entropy maxima (or minima of free en- 
ergy) are spherically symmetric. Introducing the cen- 
tral density po = p(0), the central potential $o = $(0), 
the new field ifi = /3m(<I> — $o) and the scaled distance 
£ = {SdGPmpoY^r, the Boltzmann distribution (|44|) 
can be rewritten 



P = Poe 



-V>(f) 



(46) 



Substituting this relation in the modified Poisson equa- 
tion (|43[) . we obtain the modified Emden equation 



1 d 



d-i# 



A. 



(47) 



where A = ~p/ po plays the role of the inverse central den- 
sity. Since $'(0) = for a spherically symmetric system, 
the boundary conditions at the origin are 



^(0) = V>'(0) = 0. 



(48) 



The ordinary Emden equation [53[ is recovered for A = 0, 
i.e. for very large central densities with respect to the 
average density. The function e~^^> is plotted in Figs. 
[T]and [2] for different values of A and different dimensions 
of space d. It presents an infinity of oscillations. For 
d = 1, the oscillations are undamped and their period 
is given by Eq. (|Elip . For d > 2, the oscillations are 
damped and the function ip(£.) tends to the asymptotic 
value — In A for £ — > +oo. 

We assume that the system is enclosed in a spheri- 
cal box of radius R. The normalized box radius a = 
{SdGflmpo) 1 / 2 ]! is determined by the boundary condi- 
tion $'(i?) = that becomes 



ip'(a) = 0. 



(49) 



For a given value of A, we need to integrate the modi- 
fied Emden equation (|4"T)) - (|4"5]) until the point £ = a such 
that ip'( a ) = 0- Since the function presents an in- 
finite number of oscillations, this determines an infinity 
of solutions ai(A), «2(A),... that will correspond to dif- 
ferent branches in the following diagrams. Once a n (X) 



7 




400 



5^200 




FIG. 1: The function e"^ for d = 1 and A = 0.5 < 1 (bottom) 
or A = 2 > 1 (top). In d = 1, the oscillations are undamped. 



FIG. 3: Inverse temperature r\ as a function of the central 
density 1/A for the first three branches in d = 1. 




FIG. 2: The function e - " 1 for d = 2 and A = 0.5 < 1 (bottom) 
or A = 2 > 1 (top). In d > 2, the oscillations are damped. 
The case d = 3 (not represented) is similar. 



is determined, the density profile is given by Eq. (pf()|). 
The density profile is extremum at the center and at the 
boundary. On the n-th branch, the density profile shows 
n "clusters" corresponding to the oscillations of e - '^-'. 
Close to the origin, the density increases for A > 1 while 
it decreases for A < 1. The homogeneous state ift = cor- 
responds to A = 1. This solution is degenerate because 
the boundary condition P5]l is satisfied for any a. 

Remark: When A — > 0, corresponding to large values of 
the central density, we expect to obtain results similar to 
those obtained for the usual Newtonian model since the 
differential equation (|47[) reduces to the ordinary Emden 
equation. However, the results are different because the 
boundary conditions are not the same. In the Newto- 
nian model, the force at the boundary is non zero (for 
a spherically symmetric system, according to the Gauss 
theorem, we have 3>'(i?) = GM/i? d_1 ) while in the mod- 
ified Newtonian model the force at the boundary is zero 
(<f>'(i?) = 0). Therefore, strictly speaking, the Newtonian 
and the modified Newtonian models behave differently 
even when po +oo. Nevertheless, for large central 



concentrations, the Newtonian solution provides a good 
approximation of the modified Newtonian solution in the 
core (see Appendix |E|) . 



C. The temperature 

We must now relate the normalized central density 1/A 
to the temperature T. Recalling that ~p = M /V with 
V = \S d R d , we obtain 



p dM 1 ,GMmB 1 
A = — = — — 7 — = d- 



Po S d R d po 
Introducing the normalized temperature 

f3GMm 



V 



R d ~ 



we find the relation 



V 



-Xa 



(50) 



(51) 



(52) 



Recalling that a = a n {X) for the n-th branch, this equa- 
tion gives the relation between the inverse temperature r\ 
and the central density 1/A for the n-th branch. In Figs. 
[31 S] and [SJ we plot the inverse temperature n as a func- 
tion of the central density 1/ A for the first three branches 
?i = 1,2,3 in different dimensions of space d = 1, 2, 3. 

Let us discuss the asymptotic behaviors of the tem- 
perature (we only describe the first branch n = 1) and 
compare with the Newtonian model (see, e.g., HH): 

• In d — 1: for the ordinary Newtonian model, the 
series of equilibria is parameterized by a, which is a 
measure of the central density. When a — > +oo, the 
distribution tends to a Dirac peak p = M8{x) and the 
inverse temperature rj — > +oo. When a — > 0, the dis- 
tribution is homogeneous and the inverse temperature 
r) —> 0. For the modified Newtonian model, the se- 
ries of equilibria is parameterized by the central den- 
sity 1/A. When 1/A —> +oo, the distribution tends to 
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FIG. 4: Inverse temperature r\ as a function of the central 
density 1/A for the first three branches in d = 2. 



• In d = 3: for the ordinary Newtonian model, the se- 
ries of equilibria is parameterized by a. When a — > +oo, 
the distribution tends to the singular isothermal sphere 
p s (r) = 1 / '(27T Gfimr 2 ) and the inverse temperature r\ — > 
T} s = 2. The curve 77(a) displays damped oscillations 
around this value. When a — > 0, the distribution is ho- 
mogeneous and the inverse temperature 77 — >■ 0. For the 
modified Newtonian model, the series of equilibria is pa- 
rameterized by 1/A. When 1/A — > +00, the distribution 
is concentrated at the center and we numerically find 
that 77 — > 3.05... (the value is different from the Newto- 
nian result rj s = 2 due to different boundary conditions). 
The curve 77(A) displays damped oscillations around this 
value. When A = 1, the distribution is homogeneous and 



7/ = 77* = \x\ ~ 6.7302445 (see Appendix E}. When 



1/A — > 0, the distribution is concentrated at the bound- 
ary and we numerically find that 77 — > +00. 




5 10 



15 20 



FIG. 5: Inverse temperature r\ as a function of the central 
density 1/A for the first three branches in d = 3. 



a Dirac peak p = MS(x) and 77 — > +00 with the same 
asymptotic behavior as in the Newtonian model (see Ap- 
pendix [Ej) . When A — 1, the distribution is homoge- 
neous and 77 = 77* = 7T 2 ~ 9.8696044 (see Appendix E}. 
When 1/A — s- 0, the distribution tends to a Dirac peak 
p = ^§-{5{x — R) + 5{x + R)) concentrated at the box and 



M 
2 

77 — > +00 
• In d 



2: for the ordinary Newtonian model, the 
series of equilibria is parameterized by a. When a — ¥ 
+00, the distribution tends to a Dirac peak p = M8(r) 
and the inverse temperature tends to 77 c = 4. When 
a — > 0, the distribution is homogeneous and the inverse 
temperature 77 — > 0. For the modified Newtonian model, 
the series of equilibria is parameterized by 1/A. When 
1/A — >• +00, the distribution tends to a Dirac peak p = 
M6(r) and 77 — > 77c = 4 (since the density is very much 
concentrated, the boundary conditions do not matter and 
we recover the same results as in the Newtonian case). 
When A = 1, the distribution is homogeneous and 77 = 
rf c = ij^ ~ 7.3410008 (see Appendix EJ. When 1/A -> 
0, the distribution is concentrated at the boundary and 
i] —t +00. 



D. The energy 

Wc must also relate the normalized central density 1/A 
to the energy E. The total energy is given by (see Ap- 
pendix [(J: 



E = J f V ldvdw+ l - J (p-p)$dr. 



(53) 



Using the Maxwell-Boltzmann distribution (|10p . the ki- 
netic energy is simply 



K = -Nk B T. 
2 



(54) 



Using the modified Poisson equation (|43[) and an integra- 
tion by parts, the potential energy can be written 



W 



^/(v®)>*. 



2SdG 

The total energy E = K + W is therefore given by 



E = -Nk B 
2 



(55) 



(56) 



Introducing the dimensionlcss variables defined previ- 
ously, recalling that r = £i?/a, and introducing the nor- 
malized energy 



A EE - 



ER 



d-2 



GM 2 



(57) 



we obtain 



A = 



1 1 



2i] 2rf a 



Recalling that a = a n (X) and 77 = 77™ (A) for the n-th 
branch, this equation gives the relation between the en- 
ergy A and the central density 1/A for the n-th branch. 
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FIG. 6: Energy A as a function of the central density 1/A for 
the first three branches in d — 1. 
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FIG. 8: Energy A as a function of the central density 1/A for 
the first three branches in d = 3. 
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FIG. 7: Energy A as a function of the central density 1/A for 
the first three branches in d — 2. 



In Figs. |51 [7] and [5J we plot the normalized energy 
A as a function of the central density 1/A for the first 
three branches n = 1, 2, 3 in different dimensions of space 
d= 1,2,3. 

Let us consider the asymptotic behaviors of the energy 
(we only describe the first branch n = 1) and compare 
with the Newtonian model (see, e.g., (29[): 

• In d — 1 : for the ordinary Newtonian model, the se- 
ries of equilibria is parameterized by a, which is a mea- 
sure of the central density. When a +oo, the distri- 
bution tends to a Dirac peak p = M6(x) and the energy 
A — > 0. When a — > 0, the distribution is homogeneous 
and the energy A — > — oo. For the modified Newtonian 
model, the series of equilibria is parameterized by the 
central density 1/A. When 1/A — > +oo, the distribution 

tends to a Dirac peak p = M5(x) and A — > AmLc = 1/6 
(see Appendix |D|) . When A = 1 the distribution is ho- 
mogeneous and A = A* = -1/(2??*) ~ -0.0506606. 
When 1/A — » 0, the distribution tends to a Dirac peak 
p = ^-(S(x — R) + 5(x + R)) concentrated at the box and 

A -> A (1) 

• In d = 2: for the ordinary Newtonian model, the se- 



ries of equilibria is parameterized by a. When a — > +oo, 
the distribution tends to a Dirac peak p = MS(r) and 
the energy A ^ +oo. When a — > 0, the distribution is 
homogeneous and the energy A — > — oo. For the modified 
Newtonian model, the series of equilibria is parameter- 
ized by 1/A. When 1/A — > +oo, the distribution tends 
to a Dirac peak p = MS(r) and A — » +oo. When A = 1 
the distribution is homogeneous and A = A* — —1/rj* ~ 
-0.13622121. When 1/A -> 0, the distribution is con- 
centrated at the boundary and we numerically find that 
A-> 0.1. 

• In d = 3: for the ordinary Newtonian model, the 
series of equilibria is parameterized by a. When a — > 
+oo, the distribution tends to the singular isothermal 
sphere p s (r) = 1 / (2nG (3mr 2 ) with energy A s = 1/4. The 
curve A(a) undergoes damped oscillations around this 
value. When a — > 0, the distribution is homogeneous 
and the energy A — > — oo. For the modified Newtonian 
model, the series of equilibria is parameterized by 1/A. 
When 1/A —> +oo, the distribution is concentrated at the 
center and we numerically find that A — > —0.38... (the 
value is different from the Newtonian result A s = 1/4 
due to different boundary conditions). The curve A(A) 
undergoes damped oscillations around this value. When 
A = 1 the distribution is homogeneous and A = A* = 
—3/(277*) - -0.22287452. When A -> 0, the distribution 
is concentrated at the boundary and we numerically find 
that A -> 0.05. 



E. The entropy and the free energy 

Finally, we relate the central density 1 /A to the entropy 
S and to the free energy F. Using Eqs. (@}, (JTUJ and (fTTI) . 
the entropy is given by 

S = ^Nks In T - k B [ f- In -- dr. (59) 
2 J m m 
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FIG. 9: Free energy Nk as a function of the inverse temper- 
ature rj in d = 1. Note that the branches A < 1 and A > 1 
coincide. 
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FIG. 10: Free energy Nk as a function of the inverse tem- 
perature rj in d = 2. 



Substituting Eq. ([46]) in Eq. ([59)1. and introducing the 
dimcnsionless variables defined previously, we get 



5 



Nk f 



Nm V a 



In j3 — In po 



(60) 



up to some unimportant constants. Using a = 
express po m terms of a and introduc- 
ing the normalized temperature (|51[) , we finally obtain 



S 



d-2 



Nk B 
1 1 



In rj — 2 In a 



rj a d 2 j 



(61) 



up to some unimportant constants. Using the previous 
results, this expression relates the entropy S/Nks to the 
central density 1/A. The free energy is F = E — TS. In 
the following, it will be more convenient to work in terms 
of the Massieu function J = 5 — ksfiE (by an abuse of 
language, we shall often refer to J as the free energy). 
We have 



J 



S 



Nkf 



Nkf 



77 A. 



(62) 



Using the previous results, this expression relates the free 
energy J/Nk B to the central density 1/A. 

In Figs. HI [TO] and [TTJ we have plotted the free en- 
ergy J /Nks as a function of the inverse temperature rj 
(parameterized by the central density 1/A) in d = 1, 2, 3. 
In Figs. [12j [131 [14] and HH we have plotted the entropy 
S/Nks as a function of the energy A (parameterized by 
the central density 1/A) in d = 1,2, 3. In these figures, 
the solid lines without label refer to the homogeneous 
phase. The solid lines with label n = 1 refer to the 
first inhomogencous branch. The dashed lines with label 
n = 2 refer to the second inhomogencous branch. These 





1 


1 


\ -3.15 


J. 1 1 1 1 1 1 \ 1 1 1 1 1 - 




\ 32 






- \ .—,-3.25 






AH> 0\ - 3 3 
\^ A.-1 -3.35 


7 i , i , i , i , i Nj\~ 






2.7 2.8 2.y 3 3.1 3.2 
Tl 

n= 


1 


- d=3 


i 


n=2 

i 



10 15 20 25 30 35 
1 



FIG. 11: Free energy N J k as a function of the inverse tem- 
perature rj in d = 3. 



curves will be helpful in the next section to analyze the 
phase transitions in the canonical and microcanonical en- 
sembles respectively. 

Remark: Since SS = ksPSE, the extrema of entropy 
5(A) and energy E(X) coincide. Since the series of equi- 
libria -E'(A) exhibits damped oscillations for 1/A —> +00 
in d = 3 (see Fig. [5]), this implies that the curve 5(A) 
will also exhibit damped oscillations at the same loca- 
tions. Correspondingly, S(E) will present some "spikes" 
for 1/A — > +00 in d = 3 (see inset of Fig. [To]) . Similarly, 
since 5 J = —EksSp, the extrema of free energy J(A) 
and temperature /3(A) coincide. Since the series of equi- 
libria /3(A) undergoes damped oscillations for 1/A —> +00 
in d = 3 (see Fig. [5J, this implies that the curve J(A) 
will also exhibit damped oscillations at the same loca- 
tion, and that the curve J(/3) will present some "spikes" 
for 1/A — > +00 in d = 3 (see inset of Fig. [Tl"j) . In ad- 
dition, the curve J(/3) presents a minimum for 77 ~ 24.7 
corresponding to E = 0. Similar behaviors were previ- 
ously observed in the model of sclf-gravitating fcrmions 

Mm- 
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FIG. 12: Entropy ^ as a function of energy A in d = 1. 
Note that the branches A < 1 and A > 1 coincide. 
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FIG. 13: Entropy N f. as a function of energy A in d = 2. 
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FIG. 14: Enlargement of Fig. [13] The entropies of the ho- 
mogeneous phase and inhomogeneous phase become equal at 
A = A t ~ —0.146. This corresponds to a first-order phase 
transition in the microcanonical ensemble marked by the dis- 
continuity of the slope S'(E) = 1/T. 
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FIG. 15: Entropy N ^ as a function of energy A in d = 3. 



F. Caloric curves and phase transitions 

We shall now determine the caloric curve /3(E) corre- 
sponding to the modified Newtonian model. First of all, 
we note that for the homogeneous phase, the potential 
energy W = so that the energy reduces to the kinetic 
energy. Therefore, the series of equilibria of the homoge- 
neous phase is simply 



On the other hand, eliminating A between rj n (X) and 
A„(A) given by Eqs. ([52]) and (|55|) . we get the series 
of equilibria rj n (A) for the n-th inhomogeneous branch. 
The series of equilibria contain all the critical points of 
the optimization problems ([3]) and (fl"7|) . The series of 
equilibria are the same in the canonical and microcanon- 
ical ensembles because the critical points are the same. 
They contain fully stable states (global maxima of S or 
J), metastable states (local maxima of S or J) and un- 
stable states (saddle points of S or J). The stable parts 
of the series of equilibria form the caloric curves in the 
canonical and microcanonical ensembles. We shall dis- 
tinguish the strict caloric curves formed by fully stable 
states and the physical caloric curves containing fully sta- 
ble and metastable states [66| . Metastable states are im- 
portant because they can be long-lived in systems with 
long-range interaction [5ol l5o| . The caloric curves may 
differ in CE and MCE in case of ensembles inequiva- 
lence. They arc described below in different dimensions 
of space. 

Remark: In order to determine the stable branch, we 
shall compare the entropy (in MCE) or the free energy 
(in CE) of the different solutions in competition (with 
the same values of energy or temperature). However, 
this is not sufficient because a distribution could have a 
high entropy and be an unstable saddle point. A more 
rigorous study should therefore investigate the sign of the 
second order variations of entropy or free energy for each 
critical point. But this is a difficult task that is left for 
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future works. In order to find the stable states, we shall 
use physical considerations and exploit results obtained 
in related studies. 



1. The dimension d — 1 

In Fig. [IB] we plot the series of equilibria in d = 1. 

Let us first describe the canonical ensemble (CE). The 
control parameter is the inverse temperature rj and the 
stable states are maxima of free energy J at fixed mass 
AI. The homogeneous phase exists for any value of r). It 
is fully stable for rj < rf c and unstable for r\ > rj* (see 
Sec. \W} . The first branch n = 1 of inhomogeneous states 
exists only for r/ > r/* . It has a higher free energy J than 
the homogeneous phase (see Fig. [9]) and it is fully stable. 
Secondary branches of inhomogeneous states appear for 
smaller values of the temperature but they have smaller 
values of free energy J (see Fig. ^ and they are unstable 
(saddle points of free energy). Therefore, the canonical 
caloric curve displays a second order phase transition be- 
tween homogeneous and inhomogeneous states marked 
by the discontinuity of ?M at /3 = j3*. For the inhomo- 
geneous states, there exists two solutions with the same 
temperature and the same free energy but with differ- 
ent density profiles corresponding to Ai < 1 and A2 > 1 
(see Fig. [17]). Thus, the inhomogeneous branch is de- 
generate. These two states can be distinguished by their 
central density 1/A. In conclusion: (i) for r\ < rj*, there 
is only one stable state A = 1 (homogeneous); (ii) for 
77 > 77*, there are two stable states Ai < 1 and A2 > 1 
(inhomogeneous) with the same free energy and one un- 
stable state A = 1 (homogeneous). Therefore, the central 
density 1/A plays the role of an order parameter (see 
Fig. [3]). In d = 1, there exists a fully stable equilibrium 
state for any temperature. This is consistent with the 
usual Newtonian model in d = 1 [2(], Ull ■ This is also 
consistent with results of chemotaxis since it has been 
rigorously proven that the Keller-Segel model does not 
blow up in d = 1 Q . 

Let us now describe the microcanonical ensemble 
(MCE). The control parameter is the energy A and the 
stable states are maxima of entropy S at fixed mass M 
and energy E. The homogeneous phase exists for any 
value of energy A < 0. It is fully stable for A < A* and 
unstable for A > A* (see Sec. E$. The first branch n = 1 
of inhomogeneous states exists only for A* < A < A max . 
It has a higher entropy S than the homogeneous phase 
(see Fig. IT!?]) and it is fully stable. Secondary branches 
of inhomogeneous states appear for smaller values of the 
energy but they have smaller values of entropy S (see 
Fig. IT2"j) and they are unstable (saddle points of en- 
tropy). Therefore, the microcanonical caloric curve dis- 
plays a second order phase transition marked by the dis- 
continuity of ^ at E = E*. For the inhomogeneous 
states, there exists two solutions with the same energy 
and the same entropy but with different density profiles 
corresponding to Ai < 1 and A2 > 1. Thus, the inhomo- 
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FIG. 16: Series of equilibria in d = 1. The caloric curve dis- 
plays a second order phase transition in CE and MCE taking 
place at n = rf c and A = A* (corresponding to A = 1). It is 
marked by the discontinuity of df3/dE in MCE or dE/d/3 in 
CE. Note that the branches A < 1 and A > 1 coincide. The 
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FIG. 17: Density profiles of the two stable inhomogeneous 
solutions Ai = 0.69 < 1 and A2 = 1.54 > 1 corresponding 
to X] = 10 in d = 1. We have also represented the unstable 
homogeneous solution. 



geneous branch is degenerate. These two states can be 
distinguished by their central density 1/A. In conclusion: 
(i) for A < A* , there is only one stable state A = 1 (homo- 
geneous); (ii) for A* < A < 0, there are two stable states 
Ai < 1 and A2 > 1 (inhomogeneous) with the same en- 
tropy and one unstable state A = 1 (homogeneous), (iii) 
for < A < A max , there are two stable states Ai < 1 and 
A2 > 1 (inhomogeneous) with the same entropy. There- 
fore, the central density 1/A plays the role of an order 
parameter (see Fig. [6j. In d = 1, there exists a fully 
stable equilibrium state for any accessible energy. This 
is consistent with the usual Newtonian model in d = 1 

The caloric curve, corresponding to the fully stable 
states in the series of equilibria, is denoted by (S) in 
Fig. [TBI The branch (U) corresponds to unstable states. 
There exists a fully stable equilibrium state for any ac- 
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ccssiblc values of energy in MCE and temperature in CE. 
The microcanonical and canonical ensembles are equiva- 
lent (like in the Newtonian case). 

In conclusion, the system displays second order phase 
transition in CE and MCE. This is similar to the HMF 
model dHHI. 



2. The dimension d — 2 

In Fig. [18] we plot the series of equilibria in d = 2. 

Let us first describe the canonical ensemble (CE). The 
control parameter is the inverse temperature rj. The ho- 
mogeneous phase exists for any 77. It is stable for 77 < -q* 
and unstable for r\ > rf c (see Sec. |VJ). The first branch 
n = 1 of inhomogeneous states exists for ?; > r\ c = 4 and 
it connects the homogeneous branch at r/*. For rj < rj*, 
it has a lower free energy J than the homogeneous phase 
(see Fig. ITU)) and it is unstable. For 77 > 77*, it has a 
higher free energy J than the homogeneous phase (see 
Fig. [10]). However, it is expected to be unstable or, pos- 
sibly, metastable (to settle this issue we have to study 
the sign of the second order variations of free energy as 
explained above). Secondary inhomogeneous branches 
appear for smaller values of the temperature but they 
have smaller values of the free energy (see Fig. [TO]) and 
they are unstable. The homogeneous branch is expected 
to be fully stable for rj < r\ c = 4 and metastable for 
rjc = 4 < 7] < 77* (see Fig. [19]). These conclusions 
are motivated by two arguments: (i) in the Newtonian 
model in d = 2, we know that there exists a fully sta- 
ble equilibrium state for r) < r) c = 4 and no equilibrium 
state for 77 > r\ c = 4. In that case, the system under- 
goes an isothermal collapse [IH, [2!|. For 77 > r? c = 4, 
there is no global maximum of free energy J because we 
can make it diverge by creating a Dirac peak containing 
all the particles. In the modified Newtonian model, the 
same argument applies since it is independent of bound- 
ary conditions. Since we know that the homogeneous 
branch is stable for 77 < 77*, we conclude that it must be 
metastable in the range r) c = 4 < 77 < 77* . There is there- 
fore a zeroth order phase transition at ry c = 4 marked by 
the discontinuity of the free energy, (ii) In the chemo- 
tactic literature, it has been rigorously established that 
the Keller-Segel model in d = 2 does not blow up for 
V < Vc = 4 while it can blow up for 77 > r\ c = 4 Q. This 
is consistent with our stability results. 

Let us now describe the microcanonical ensemble 
(MCE). The control parameter is the energy A. The 
homogeneous phase exists for all values of A < 0. It is 
stable for A < A* and unstable for A > A* (see Sec. 
|V]) . The first branch n = 1 of inhomogeneous states ex- 
ists for A > A* and it connects the homogeneous branch 
at A*. We see that the inhomogeneous branch (3(E) is 
multi- valued. Considering the value of the entropy in the 
different phases (see Figs. [T3l and IT4"]) . the caloric curve 
is expected to display a microcanonical first order phase 
transition at A = A f ~ —0.146 marked by the discontinu- 
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FIG. 18: Series of equilibria in d = 2. The first inhomoge- 
neous branch n = 1 tends to a plateau rj c = 4 for large central 
densities 1/A — > +00 due to the formation of a Dirac peak. 
This is similar to the plateau appearing in the caloric curve 
of the classical self-gravitating gas [2^ . 
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FIG. 19: Caloric curve in the canonical ensemble in d — 2. 
The homogeneous branch is fully stable for ?7 < rjc — 4. 
metastable for rj c < n < rj* and unstable for 77 > 77* . The inho- 
mogeneous branch is always unstable (or, possibly, metastable 
for 77 > 77*). For sufficiently low temperatures, the system can 
experience an isothermal collapse. 



ity of the temperature (see Fig. [20]). The energy of tran- 
sition has been determined by comparing the entropy of 
the homogeneous and inhomogeneous phases and looking 
at which point the curves S(E) intersect (see Fig. IT4"]) . 
Equivalently, it can be obtained by performing a vertical 
Maxwell construction [3l[. The homogeneous phase is 
fully stable for A < A t , metastable for A t < A < A* and 
unstable for A > A*. The lower part of the first inhomo- 
geneous branch is fully stable for A > A t and metastable 
for A* < A < A t . The upper part of the first inhomoge- 
neous branch is unstable for A* < A < A* . For A > A* , it 
is unstable or, possibly, metastable. Secondary inhomo- 
geneous branches appear for smaller values of the energy 
but they have smaller values of the entropy (see Fig. [T5]) 
and they are unstable. The stable states of the inhomo- 
geneous branch have 1/A > 1 indicating that the density 
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FIG. 20: Caloric curve in the microcanonical ensemble in 
d — 2. A first-order phase transition is expected to take 
place in the microcanonical ensemble at A = At. Note that 
the lower branch has negative specific heats. A* and possibly 
A* represent microcanonical spinodal points marking the end 
of the metastable phase. 



is concentrated at the center. The possibly metastable 
states for A > A* have 1/A < 1 indicating that the den- 
sity is concentrated near the box. In conclusion, there 
exists a fully stable equilibrium state for any value of en- 
ergy. This is similar to the Newtonian model in d = 2 
24, [2^]. However, in the modified Newtonian model, we 
expect a first order phase transition at A t that is not 
present in the Newtonian model. 

The strict caloric curve, corresponding to the fully sta- 
ble states (global maxima) in the series of equilibria, is 
denoted (S) in Figs. [T9land[20l The unstable states (sad- 
dle points) arc denoted (U) and the metastable states 
(local maxima) are denoted (M) . There exists a fully sta- 
ble equilibrium state for any accessible value of energy in 
MCE and for sufficiently high values of the temperature 
in CE (77 < f]c = 4). Here, the microcanonical and canon- 
ical ensembles are inequivalent (unlike in the Newtonian 
case). In particular, the lower part of the first inhomo- 
geneous branch is stable in MCE while it is unstable in 
CE. This branch has negative specific heats C < (see 
Fig. [2U)| which is not possible in the canonical ensemble. 

In conclusion, the system displays a zeroth order phase 
transition in CE (associated with an isothermal collapse) 
and a first order phase transition in MCE. Note also that 
the energy E{0) and its first derivative E'((3) are contin- 
uous at the critical point /3* but its second derivative 
E"(/3) is discontinuous. Provided that the inhomoge- 
neous branch for 77 > 77* is metastable, this would corre- 
spond to a third order canonical phase transition between 
a homogeneous metastable state and an inhomogeneous 
metastable state. 



Let us first describe the canonical ensemble (CE). The 
control parameter is 77. The homogeneous phase exists 
for all 77. It is stable for 77 < 77* and unstable for 77 > 77* 
(see Sec. |V|) . The first branch n = 1 of inhomogeneous 
states exists for 77 > 2.64 and it connects the homoge- 
neous branch at 77 = 77*. For large central densities 1/A, 
it forms a spiral towards a singular solution. For 77 < 77* , 
it has a lower free energy J than the homogeneous phase 
(see Fig. [TT|) and it is unstable. For 77 > 77*, it has a 
higher free energy J than the homogeneous phase (see 
Fig. [TT]). However, it is expected to be unstable or, pos- 
sibly, metastable. Secondary inhomogeneous branches 
appear for smaller values of the temperature but they 
have a higher value of free energy J and they arc unsta- 
ble. The homogeneous branch is metastable for 77 < 77*. 
These conclusions are motivated by two arguments: (i) 
in the Newtonian model in d = 3, we know that there is 
no fully stable equilibrium state in CE. The system can 
undergo an isothermal collapse [29j . There is no global 
maximum of free energy J because we can make it diverge 
by creating a Dirac peak containing all the particles [2l| . 
In the modified Newtonian model, the same argument 
applies since it is independent on boundary conditions. 
Since we know that the homogeneous branch is stable 
for 77 < 77*, then it can only be metastable. (ii) In the 
chemotactic literature, it has been rigorously established 
that the Keller-Segel model in d = 3 can blow up for any 
77 Q. This is consistent with our stability results. 
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3. The dimension d = 3 



In Fig. [21] we plot the series of equilibria in d = 3. 



FIG. 21: Series of equilibria in d = 3. The inhomogeneous 
branch forms a spiral for large central densities 1/A — » +00 
due to the damped oscillations of the inverse temperature 77(A) 
and energy A(A). This is similar to the spiral appearing in 
the series of equilibria of the classical self-gravitating gas as 
we approach the singular isothermal sphere [29| . 



Let us now describe the microcanonical ensemble 
(MCE). The control parameter is the energy A. The 
homogeneous phase exists for all A < 0. It is stable for 
A < A* and unstable for A > A* (see Sec. EJ. The 
first branch re = 1 of inhomogeneous states exists for 
A > —0.405 and it connects the homogeneous branch at 
A = A*. For large central densities 1/A, it forms a spiral 
towards a singular solution. For A < A*, it has a lower 
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entropy S than the homogeneous phase (see Fig. [T5|) 
and it is unstable. For A > A*, it has a higher entropy 
than the homogeneous phase (see Fig. [T5|) . However, it 
is expected to be unstable or, possibly, metastable. Sec- 
ondary inhomogeneous branches appear for smaller val- 
ues of the energy but they have a lower value of entropy 
S and they are unstable. The homogeneous branch is 
metastable for A < A*. These conclusions are motivated 
by two arguments: (i) in the Newtonian model in d = 3, 
we know that there is no fully stable equilibrium state in 
MCE. The system undergoes a gravothermal catastrophe 
[l3l [f| . There is no global maximum of entropy S at 
fixed mass and energy because we can make it di verge by 
creating a binary star surrounded by a hot halo [22|, [3l| . 
In the modified Newtonian model, the same argument 
applies. Since we know that the homogeneous branch is 
stable for A < A*, then it can only be metastable. 

There is no strict caloric curve since there is no fully 
stable states (global maxima). But there is a physical 
caloric curve made of metastable states (local maxima) 
denoted (M) in Fig. [2TJ The unstable states (saddle 
points) are denoted (U). Here, the microcanonical and 
canonical ensembles, regarding the metastable states, are 
equivalent unlike in the Newtonian case. This is because 
the homogeneous branch and the inhomogeneous branch 
connect each other at a single point at A = 1 by making 
a cusp (see inset in Fig. l2"Tj) while the Newtonian series 
of equilibria is smooth and presents two distinct turning 
points of temperature and energy (denoted CE and MCE 
in Fig. 8 of [3l| ) separated by a region of negative specific 
heats. 

In conclusion, if we take metastable states into ac- 
count, the system displays a zeroth order phase tran- 
sition in CE and MCE corresponding to a discontinuity 
of entropy or free energy. They are associated with an 
isothermal collapse or a gravothermal catastrophe respec- 
tively. 

Remark: There is no natural external parameter in the 
modified Newtonian model. However, the dimension of 
space d could play the role of an effective external pa- 
rameter. The preceding results predict the existence of a 
critical dimension d c between 1 and 2 at which the phase 
transition passes from second order (d < d c ) to first order 
(d > d c ). However, this transition turns out to occur in 
a very small range of parameters since we find that the 
critical dimension d c is between 1 and d = 1.00001 and 
the concerned range of energies and temperatures is ex- 
tremely narrow. We have not investigated this transition 
in detail since the dimension of space is not a physical 
(tunable) parameter. Furthermore, in the next model, we 
have an external parameter /1 played by screening length 
that is more relevant. 



IV. THE SCREENED NEWTONIAN MODEL 

In this section, we discuss phase transitions that ap- 
pear in the screened Newtonian model corresponding to 



an attractive Yukawa potential. 

A. Physical motivation of the model 

We consider a system of particles interacting via the 
potential $(r,i) that is solution of the screened Poisson 
equation 



kl<S> 



SdGp, 



(64) 



where ko is the inverse of the screening length. At statis- 
tical equilibrium, the density is given by the Boltzmann 
distribution 



p = Ae~ f3m *. 



(65) 



We assume that the system is confined in a finite domain 
(box) and we impose the Neumann boundary conditions 



V$n = 0, Vp-n = 0, 



(66) 



where n is a unit vector normal to the boundary of the 
box (the explicit expression of the potential in d = 1 is 
given in Appendix [Bjl . This model admits spatially ho- 
mogeneous solutions (p = po and $ = $ with — k^ n = 
SdGpo) at any temperature. It also admits spatially inho- 
mogeneous solutions at sufficiently low temperatures. We 
shall study this model in arbitrary dimensions of space 
d with explicit computations for d = 1,2,3. This model 
has different physical applications: 

(i) It describes a system of particles interacting via a 
screened attractive (Newtonian) potential. 

(ii) By a proper reinterpretation of the parameters, 
the field equation (|64| describes the relation between the 
concentration of the chemical and the density of bacteria 
in the Keller-Segel model (|57|) where the degradation of 
the chemical reduces the range of the interaction. In 
that case, the boundary conditions are of the form (|66|) . 
Furthermore, the relevant ensemble is the CE since the 
KS model has a canonical structure. This model has been 
studied by Childress & Percus [HH in d = 1 using an 
approach different from the one we are going to develop. 

For the sake of generality, we shall study this model 
in the microcanonical and canonical ensembles in any 
dimension of space. 



B. The screened Emden equation 

In the screened Newtonian model, the equilibrium den- 
sity profile is given by the Boltzmann distribution (|65j) 
coupled to the screened Poisson equation As in 

Sec. IIIIB1 we look for spherically symmetric solutions. 
Introducing the central density po = p(0), the central po- 
tential $0 = $(0), the new field tp = Pm(^ - $ ) and 
the scaled distance £ = (SdG/Smpo^^r the Boltzmann 
distribution can be rewritten 



-V>(«) 



(67) 
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FIG. 22: Inverse temperature rj as a function of a for the first 
two branches in d = 1. We have taken fi = 1 < fj, c . 
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FIG. 23: Inverse temperature rj as a function of a for the first 
two branches in d = 1. We have taken fx = 10 > /i c - 



Substituting this relation in the screened Poisson equa- 
tion (|64|). we obtain the screened Emden equation 

where k = k a / (SdGfimpo) 1 ^ 2 and A = —kQ^> /SdGp - 
The boundary conditions at the origin are 



7/Ho) = vAo) = o. 



(69) 



The normalized box radius is a — (SdGfimpo) 1 / 2 ]! and 
the boundary condition $'(i?) = becomes 

ip'(a) = 0. (70) 

Introducing the normalized screening length 

fi = k R, (71) 

the parameter k can be rewritten n = p/a. For given /z, 
we solve the problem as follows: (i) We fix a. (ii) K = 
/i/a is then given, (iii) We determine A by an iterative 
method such that ip'( a ) — 0- (i y ) We obtain different 
solutions A„(q;) determining different branches n = 1, 
n = 2 etc. This procedure determines for each value of 
a, and for each branch, the normalized density profile 
g— The homogeneous solution corresponds to ip = 
and A = 1. This solution is degenerate because the 
boundary condition (|70p is satisfied for any value of a. 



C. The temperature 

We must now relate the parameter a to the tempera- 
ture T . Introducing the dimcnsionlcss variables defined 
previously and recalling that r = Rt;/a, the mass can be 
written 
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FIG. 24: Inverse temperature rj as a function of a for the first 
two branches in d — 2. For a — > +oo, the inverse temperature 
of the first branch tends to r\ c = 4. We have taken fi = 1. 




FIG. 25: Inverse temperature as a function of a for the first 
two branches in d = 3. For a — > +oo, the inverse temperature 
of the first branch undergoes damped oscillations around the 
value rj s ~ 3.25. We have taken fj, = 1. 
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FIG. 26: Energy A as a function of a for the first two branches 
in d = 1. We have taken fi = 1 < /i c . 



0.03 
0.025 - 
0.02 
< 0.015 
0.01 
0.005 

°0 




50 



100 150 

a 



2< )() 



250 



FIG. 28: Energy A as a function of a for the first two branches 
in d = 1. We have taken ^ = 15 > fi m - 
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FIG. 27: Energy A as a function of a for the first two branches 
in d = 1. We have taken p = 10 > /i c . 




FIG. 29: Energy A as a function of a for the first two branches 
in d = 2. We have taken p = 1. 



D. The energy 



Using a — {SdCfirap^Y^R and introducing the dimen- 
sionlcss temperature (|51[) . we obtain 



^^[^•^dd- (73) 



This equation gives the relation between the inverse tem- 
perature r\ and a for the n-th branch. In Figs. 1221 1231 
1241 and 1251 we plot the inverse temperature rj as a func- 
tion of a for the first two branches n = 1, 2 in different 
dimensions of space <f = 1, 2, 3. The discussion is similar 
to the one given in Sec. IIII CI We have also represented 
the branch corresponding to the homogeneous solution. 
Its equation is given by r\ = a 2 /d. The branch n = 1 of 
inhomogeneous solutions connects the branch of homo- 
geneous solutions at a 2 c = drj* (see Appendix IF)) . 



We must also relate a to the energy E. The total 
energy is given by 

E = J f— drdv + \J dr. (74) 

Using the Maxwcll-Boltzmann distribution (|10|) . the ki- 
netic energy is simply 

K = ^Nk B T. (75) 

Using the screened Poisson equation (|64p and integrating 
by parts, the potential energy can be written 

W - -— ^ J [(V$) 2 + fc 2 $ 2 ] dr. (76) 

The total energy E = K + W is therefore given by 

E = ^Nk B T - _L J [(V$) 2 + fc 2 * 2 ] dr. (77) 
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FIG. 30: Energy A as a function of a for the first two branches 
in d = 3. For a — > +oo, the energy of the first branch un- 
dergoes damped oscillations around the value A s ~ 1.13. We 
have taken fj, = 1. 
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FIG. 31: Free energy N ^ as a function of the inverse tem- 
perature r\ for d = 1. We have taken fi = 1 < /i c . 



Introducing the dimensionlcss variables defined previ- 
ously, recalling that r = (,R/a and fi = k$R, and in- 
troducing the normalized energy (|57[) . we obtain 

2 7? + 2r ? 2 a d - 2 J \d£ J * 4 

T^ + ^o) 2 ^- 1 ^. (78) 

Using the expressions of k and A following Eq. (|68p . we 
find that 

/3m$ = -4- ( 79 ) 

so that, finally, 

This equation gives the relation between the energy A 
and a for the rc-th branch. In Figs. |2"rJl 1271 1251 and 
[501 we plot the energy A as a function of a for the first 
two branches n = 1 and n = 2 in different dimensions 
of space d = 1,2,3. The discussion is similar to the 
one given in Sec. MI Dl We have also represented the 
branch corresponding to the homogeneous solution. Us- 
ing Eq. ([87]) and r\ = a 2 /d, its equation is given by 
A = -d 2 /(2a 2 ) +d/(2fi 2 ). 

E. The entropy and the free energy 

Finally, we relate a to the entropy S and to the free 
energy F. The entropy is given by 

S= ~Nk B lnT-k B i — \n — dx. (81) 
2 J m m 
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FIG. 32: Free energy N J k as a function of the inverse tem- 
perature rj for d = 1. We have taken fj, = 10 > /i c . The 
free energies of the homogeneous phase and inhomogeneous 
phase become equal at r\ — r\ t fju) . This corresponds to a first 
order phase transition in the canonical ensemble marked by 
the discontinuity of the slope J'(/3) = — E. 
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FIG. 33: Free energy Nk as a function of the inverse tem- 
perature r\ for d — 2. We have taken fj, = 1. 
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FIG. 34: Free energy N J k as a function of the inverse tem- 
perature rj for d — 3. We have taken p = 1. 



FIG. 37: Entropy N j_ as a function of the energy A for d = 2. 
We have taken /i = 1. 
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FIG. 35: Entropy as a function of the energy A for d = 1. 
We have taken p = 1 < p c . 
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FIG. 38: Entropy N ^ as a function of the energy A for d = 3. 
We have taken p = 1. 
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FIG. 36: Entropy jJr as a function of the energy A for d = 1. 
We have taken p — 10 > /i c . There is a (small) convex dip 
associated with the region of negative specific heats in the 
microcanonical ensemble. 



We can proceed exactly as in Sec. IIII El and obtain 
S d-2 



Nk B 
1 1 



In T) — 2 In a 



Ve"'''^- 1 d£, 



(82) 



up to unimportant constants. However, we can 
also obtain a simpler expression. Substituting p = 
poe -/3m(*-<&o) in Eq (JsiJ)^ we obtain 

S =-Nk B \nT-k B [ f-ln^-dr 
2 mm 



+k B f3 J p($-$ )rfr- 



This can be rewritten 
S d 



Nk B 



9 R F 1 

In - In po + — - /3m$ , 



(83) 



(84) 



up to unimportant constants. Finally, using Eqs. 
(|T9"|) . (|FT|) . (|57|) and the relations n — p/a and a — 
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(SdG(3mp ) 1 ^ 2 R, we obtain 



d-2 



Nk B 



In 77 — 2 In a — 2 A77 



\a 2 



(85) 



which docs not involve new integrals. Using the previous 
results, this expression relates the entropy S/Nks to a. 
The free energy is F = E — TS. In the following, it 
will be more convenient to work in terms of the Massicu 
function J = S — ksPE (by an abuse of language, we 
shall often refer to J as the free energy). We have 



J 



S 



Nk f 



Nk f 



T]A. 



(86) 



Using the previous results, this expression relates the free 
energy J/Nk B to a. 

In Figs. El El ESI and [Ml we have plotted the free 
energy J/NkB as a function of the inverse temperature 
77 (parameterized by a) in d = 1,2,3. In Figs. ESI EH 
l37l and EE]), we have plotted the entropy S/Nks as a 
function of the energy A (parameterized by a) in d = 
1,2,3. 



F. Caloric curve 

We shall now determine the caloric curve (3(E) corre- 
sponding to the screened Newtonian model. First of all, 
we note that, for the homogeneous phase, we have p = p a 
and $ = $0 with — k 2 $ n = SdGpo (or cquivalently V> = 0, 
A = 1 and a 2 = drf). Therefore, the relationship between 
the energy and the temperature can be written 



A 



d d 
2^ + V' 



(87) 



A, 



d/(2p 2 



This shows that 77 — > +00 for A 
On the other hand, eliminating a between rj n (a) and 
A n (a) given by Eqs. (|T3"|) and we get the series 

of equilibria for the 77,-th inhomogencous branch. The 
series of equilibria (critical points) and the caloric curves 
(stable states) in CE and MCE arc described below for 
different dimensions of space. 



1. The dimension d = 1 

In Figs. ESI and 001 we plot the series of equilibria in 
d = 1 for different values of the screening parameter p. 

Let us first describe the canonical ensemble (CE). The 
control parameter is the inverse temperature 77 and the 
stable states are maxima of free energy J at fixed mass 
M. The homogeneous phase exists for any value of 77. 
It is stable for 77 < 77* and unstable for 77 > 77* (see 
Sec. EJ. Comparing Figs. ESI and FUJI we see that the 
screened Newtonian model is characterized by a pitchfork 
bifurcation at 77 = rj* . The pitchfork bifurcation is super- 
critical if p < p c = \/2~7T ~ 4.4428829 and sub-critical if 



p > p c . This interesting transition was first evidenced by 
Childress & Percus [HI using a different approach. In our 
thermodynamical approach, this implies the existence of 
a canonical tricritical point at p c = ^/2tt. For p < p c 
the phase transition is second order and for p > p c the 
phase transition is first order. 

Let us first consider p < p c (see Fig. ESJ- The discus- 
sion is similar to that given for the modified Newtonian 
model. The first branch n = 1 of inhomogencous states 
exists only for 77 > 77* . It has a higher free energy J than 
the homogeneous phase (see Fig. EH) and it is fully sta- 
ble. Secondary branches appear for smaller values of the 
temperature but they have smaller values of free energy J 
(see Fig. EH) and they are unstable (saddle points of free 
energy). Therefore, the canonical caloric curve displays 
a second order phase transition between homogeneous 
and inhomogencous states marked by the discontinuity 
of ijM at /3 = /3*. We note that, for the inhomogeneous 
states, there exists two solutions with the same tempera- 
ture and the same free energy but with different density 
profiles corresponding to a\ < a c and ct2 > ol c , where 
a c is the value of a at the point of contact with in the 
homogeneous branch. Thus, the inhomogencous branch 
is degenerate. These two states can be distinguished by 
their central density a. Since ~p/po = drj/a 2 , the solution 
Qi < a c corresponds to po < ~p and the solution ct2 > a c 
corresponds to po > ~p. The density profiles are similar 
to those represented in Fig. [17] for the modified Newto- 
nian model. In conclusion: (i) for 77 < 77*, there is only 
one stable state (homogeneous) ; (ii) for 77 > 77* , there are 
two stable states a\ < a c and a-i > a c (inhomogeneous) 
with the same free energy and one unstable state (homo- 
geneous). Therefore, the central density a plays the role 
of an order parameter (see Fig. 
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FIG. 39: Series of equilibria in d = 1 for p = 1 < p c . The 
caloric curve displays a second order phase transition in CE 
and MCE taking place at 77 = 77* and A = A* . It is marked by 
the discontinuity of dE/d/3 in CE or d/3/dE in MCE. Note 
that the branches a < a c and a > a c coincide. For p < p c , 
the CE and MCE ensembles are equivalent. 



Let us now consider p > p c (see Fig. |4U]) . The first 
branch n = 1 of inhomogencous states exists only for 
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FIG. 40: Series of equilibria in d = 1 for ji = 10 > ji c - 
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FIG. 41: Canonical caloric curve in d = 1 for /j = 10 > /i c . 
It displays a canonical first-order phase transition marked by 
the discontinuity of the energy at rj — rj t (ji). The region of 
negative specific heats is unstable in the canonical ensemble 
and replaced by a phase transition (Maxwell plateau). The 
temperatures rj* and 77* represent canonical spinodal points 
marking the end of the metastable phase. 



77 > 77* (/i). The caloric curve displays a canonical first 
order phase transition at rjt(ji) marked by the disconti- 
nuity of the energy E (see Fig. iTTj) . The temperature 
of transition 77* (yti.) can be obtained by plotting the free 
energy of the two phases as a function of the tempera- 
ture and determining at which temperature they become 
equal (see Fig. [32]). Equivalently, it can be obtained by 
performing a horizontal Maxwell construction [31| . The 
homogeneous phase is fully stable for 77 < r/ t , metastable 
for 774 < 77 < 77* and unstable for rj > rj*. The right 
branch of the inhomogeneous phase is fully stable for 
77 > rjt and metastable for 77* < 77 < rj t . The left branch 
is unstable. Note that this branch has negative specific 
heats which is not permitted in the canonical ensemble. 
Secondary branches appear for smaller values of the tem- 
perature but they have smaller values of free energy J 
and they are unstable. We also note that the branch 
of inhomogeneous states is degenerate since the curves 
a < a c and a > a c coincide. In conclusion: (i) for 
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FIG. 42: Canonical phase diagram in d = 1 exhibiting a tri- 
critical point at fi c = v27T — 4.44 and 77 ~ 29.6. We have 
represented 77c, 77* and rjt as a function of /_t. The region be- 
tween 77, and 77* contains stable and metastable states. 
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FIG. 43: Microcanonical caloric curve in d = 1 for /i = 10 > 
fXc- It displays a microcanonical second order phase transition 
marked by the discontinuity of ^ at E = E*. For ji > [i c , 
there exists a region of negative specific heats that is stable 
in the microcanonical ensemble. 



77 < 77*, there is only one stable state (homogeneous); (ii) 
for < 77 < i]* , there are three stable states (one homo- 
geneous and two inhomogeneous) and two unstable states 
(inhomogeneous); (hi) for rj > 77*, there are two stable 
states (inhomogeneous) and one unstable state (homo- 
geneous). The pairs of inhomogeneous states have the 
same free energy. Therefore, the central density a plays 
the role of an order parameter (see Fig. [23|) . 

The canonical phase diagram is represented in Fig. 02] 
where we have plotted 77*, 77, and 774 as a function of ji. 
The three temperatures coincide at the tricritical point 
fi = fi c . At that point, the phase transition goes from 
second order (/_t < fi c ) to first order (/1 > ji c ). 

The strict caloric curve (see Figs. [33] and HJ), cor- 
responding to the fully stable states, is denoted (S). 
The physical caloric curve should take into account the 
metastable states (M) because they are long-lived. The 
states (U) are unstable. We see that there exists a fully 
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FIG. 44: Microcanonical caloric curve in d = 1 for /i = 15 > 
fi m . It displays a microcanonical first order phase transition 
marked by the discontinuity of T at E = Et . The energies E* c 
and £7* are spinodal points marking the end of the metastable 
branches. Note that this first order phase transition occurs 
in an extremely small range of energies. 



FIG. 46: Microcanonical phase diagram in d = 1. We have 
represented A*, A', Ai and A2 as a function of /j,. These 
energies coincide for fi c = t/2tt ~ 4.44 and A ~ 0.0084. These 
energies delimitate respectively the region of negative specific 
heats and the region of strict ensembles inequivalence (see 
main text): the energies in these regions cannot be reached 
in the canonical ensemble. 
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FIG. 45: Microcanonical phase diagram in d = 1 exhibiting 
a tricritical point at /i m ~ 11.8 and A ~ 2.37 10 -4 . We have 
represented A c , A* and At as a function of /i. The region 
between A* and A* contains stable and metastable states. 
We again emphasize the small range of energies where this 
first order phase transition takes place. 



stable equilibrium state for any temperature and any 
screening length. This is consistent with the usual New- 
tonian model in d = 1 [2(| Hj|. This is also consistent 
with the results of chemotaxis since it has been estab- 
lished rigorously that there is no blow up in d = 1 Q . 

Let us finally describe the microcanonical ensemble 
(MCE). The control parameter is the energy A and the 
stable states arc maxima of entropy S at fixed mass M 
and energy E. The homogeneous phase exists for any 
A < A-max = d/(2fi 2 ). It is stable for A < A* and unsta- 
ble for A > A* (see Sec. EJ. Comparing Figs. H31 and 
l44l we see that there exists a microcanonical tricritical 
point at \i m ~ 11.8 and A ~ 2.37 10~ 4 (corresponding 
to ?; ~ 149.1096). For fi < fi m the phase transition is 
second order and for fi > /i m the phase transition is first 



order. 

Let us first consider /i < \i m (see Figs. [39l and l43l) . 
The first branch n — 1 of inhomogeneous states exists 
for A* < A < AmLx = l/(2//tanh(//)) (see Appendix [D]) . 
It has a higher entropy S than the homogeneous phase 
and it is stable. Secondary branches appear for smaller 
values of the energy but they have smaller values of en- 
tropy and arc unstable. The microcanonical caloric curve 
displays a second order phase transition marked by the 
discontinuity ^ at E = E*. For \i < fj, c , the specific 
heat is always positive. In that case, the microcanonical 
and canonical ensembles are equivalent. For /i > fj. c , a 
region of negative specific heats appears. This leads to 
a convex dip in the entropic curve S(E) (see Fig. [36]) . 
In that case, the microcanonical and canonical ensembles 
are inequivalent: the states with negative specific heats 
are stable in MCE while they are unstable in CE (com- 
pare Figs. ItTI and r4"3"]) . Therefore, these energies cannot 
be achieved in a canonical description. 

Let us now consider /i > fj, m (see Fig. [44)) . The first 
branch n = 1 of inhomogeneous states exists only for 
A > A*(/i). The caloric curve displays a microcanoni- 
cal first order phase transition at At(fi) marked by the 
discontinuity of the temperature T and the existence of 
metastable states. The energy of transition A t (/i) can 
be obtained by plotting the entropy of the two phases 
as a function of the energy and determining at which en- 
ergy they become equal. Equivalcntly, it can be obtained 
by performing a vertical Maxwell construction [3l1 | . The 
discussion is similar to that given in the canonical en- 
semble except that the axis are reversed. In conclusion: 
(i) for A < A*, there is only one stable state (homoge- 
neous); (ii) for A* < A < A*, there are three stable states 
(one homogeneous and two inhomogeneous) and two un- 
stable states (inhomogeneous); (iii) for A > A*, there 
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are two stable states (inhomogcncous) and one unstable 
state (homogeneous). The pairs of inhomogeneous states 
have the same entropy. Therefore, the central density a 
plays the role of an order parameter (see Fig. |2"5)) . 

The microcanonical phase diagram is represented in 
Figs. 1451 and l46l where we have plotted A*. A* and A t 
as a function of /i. The three energies coincide at the 
microcanonical tricritical point [i = /i m . At that point, 
the phase transition goes from second order (/i < /Lt m ) to 
first order (/i > /z m ). We have also represented the region 
of negative specific heats which appears at the canonical 
tricritical point /i = /Lt c . For /i c < /i < yu m , it is delimited 
by A* and A' and for /i > /i m , it is delimited by A* and 
A'. This region of negative specific heats also defines the 
physical region of ensembles inequivalence, i.e. the states 
that are stable in MCE but unstable in CE (metastable 
states are considered here as stable states). Finally, we 
have represented the strict region of ensembles inequiva- 
lence, i.e. the states that are stable in MCE but unstable 
or metastable in CE. It is delimited by Ai and A2 and, 
of course, contains the negative specific heats region. 

The strict caloric curve (see Figs. [39l l43l and [44 ]) . cor- 
responding to the fully stable states, is denoted (S). The 
states (U) are unstable. The states (M) are metastable 
but they are long-lived. We see that there exists a fully 
stable equilibrium state for any accessible energy and any 
screening length. This is consistent with the usual New- 
tonian model in d = 1 [H H<|. 

In conclusion, for p, < the system displays canon- 
ical and microcanonical second order phase transitions. 
For ii c < /i < fj, m (canonical tricritical point), the system 
displays canonical first order phase transitions and micro- 
canonical second order phase transitions. For [i > /Lt m 
(microcanonical tricritical point), the system displays 
canonical and microcanonical first order phase transi- 
tions. Note that the canonical and microcanonical tri- 
critcal points do not coincide as also observed in other 
models SI III, El. 



2. The dimensions d = 2 and d — 3 

In Figs. H7] and 051 we plot the series of equilibria in 
d = 2 and d = 3. We have considered different values of fi 
but only the case // = 1 is shown. We have observed that 
the shape of the diagrams does not significantly depend 
on the value of the screening parameter p,. Therefore, the 
description of these diagrams is similar to the one given 
in Sees. IIIIF2I and IIII F 31 for the modified Newtonian 
model. 



STABILITY OF THE HOMOGENEOUS 
PHASE 



In this section, we study the stability of the homoge- 
neous phase in the case where the potential satisfies the 
modified Poisson equation (|4"3")l or the screened Poisson 
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FIG. 47: Caloric curve in d = 2 for fj, = 1. 




FIG. 48: Caloric curve in d — 3 for fj, = 1. 



equation ([Mf. We first consider the spectral stability of 
the homogeneous phase with respect to the Smoluchowski 
equation or, equivalently, with respect to the Keller-Segel 
model. This will allow us to determine the growth rate 
(unstable case) or the damping rate (stable case) of the 
perturbation. Then, we investigate the dynamical and 
thermodynamical stability of a larger class of systems by 
determining whether the homogeneous phase is a maxi- 
mum of entropy at fixed mass and energy in MCE or a 
minimum of free energy at fixed mass in CE. 



A. Spectral stability 

We consider the mean field Smoluchowski equation 
dp 



dt 



= V • 



I (k R T. 



coupled to the modified Poisson equation (|4"3"]l or to the 
screened Poisson equation (|64)). The boundary condi- 
tions are given by Eq. (|4"5]) . Up to a change of notation, 
these equations also describe the Keller-Segel model (f3"Tj) - 
([371) or ([SI])-®- In the modified Newtonian model, the 
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homogeneous steady state satisfies 

p = p, $ = 0. (89) 

In the screened Newtonian model, it satisfies 

- fc 2 $ = s d Gp, (90) 

where p and $ are uniform. In both models, the lin- 
earized equations can be written 



£-f - -Z-A8 P + pA6$, 
at m 



A5<S> - fcg£$ = S d GSp, 



(91) 



(92) 



where ko = in the modified Newtonian model and ko ^ 
in the screened Newtonian model. 

In an infinite domain, the spectral stability of the 
homogeneous solutions of the mean field Smoluchowski 
equation (and its generalizations) coupled to Eqs. (|4"5]l 
and has been studied by Chavanis & Sire [57| who 
stressed the analogy with the Jeans instability in astro- 
physics [58[. Here, we describe how the results are mod- 
ified in a bounded domain with the boundary conditions 
P5|) . This problem was considered by Keller & Segel Q 
in d = 2. We shall perform the stability analysis in d 
dimensions. Let us call '0fe( r ) the eigenfunctions of the 
Laplacian and —k 2 the corresponding eigenvalues. They 
arc solution of 



with 



VV>fc • n = 0, 



(93) 



(94) 



on the boundary. It is easy to check that the eigenval- 
ues are necessarily negative (hence the notation —k 2 ). 
Indeed, multiplying Eq. by ifj^, integrating on the 

whole domain, and using an integration by parts, we get 
j '(V-0/c) 2 dr = k 2 J ij}\dr, which proves the result. In 
a bounded domain, their values are "quantized" (see be- 
low) . The lowest non-zero value of k will play a particular 
role as it determines the critical temperature below which 
the homogeneous phase becomes unstable. The expres- 
sion of the eigenfunctions and eigenvalues depends on the 
domain shape and on the dimension of space. In the fol- 
lowing, we shall work in a spherical box in d = 1,2,3 
dimensions. 

• In d = 1, we have 



• In d = 2, we have 

tpni = JniKir) cos(n6 

with 



kn 



tni 

IT 



(97) 



(98) 



where n is an integer and 7„i is the i-th zero of J' n (x). 
The smallest non zero eigenvalue is fcoi = loi/R where 
7oi = in = 3.83171... is the first zero of J (x) = — J\{x). 
The axisymmetric mode (n = 0) is 



ipoi = Jo(k 0l r) 
• In d = 3, we have 
1 



(99) 



Vw = -^Ji+i(hir)Yhn(9,<t>), (100) 



with 



, Hi 
kij = — 
R 



(101) 



where l,m are integers with I > \m\ and is the z-th 
zero of 



X ^'l+l/2^ X ) _ 1 _ q 



(102) 



The smallest non zero eigenvalue is fc i = 701/-R where 
7oi = x\ = 4.49341... is the first root of tan(x) = x. The 
spherically symmetric mode (I, m = 0) is 



ipo 



sin(koir) 



(103) 



The solutions of the linearized equations (|9Tj) and ([92 
can be expanded on the eigenmodes, writing 



Sp(r,t)=J2^e akt/i Mv), 



(104) 



(105) 



where the sum runs on the (quantized) eigenvalues. Sub- 
stituting Eqs. |P3gD and (fTU5|l in Eqs. $5TJ and ([52]). we 
obtain the algebraic equations 



CTfe 



^Ik 2 )A k+P k 2 B k = 0, 



(106) 



with 



ip n = cos(fc„x), 



k n — ^T7) 
K 



(95) 



(96) 



where n is an integer. The smallest non zero eigenvalue 
is k\ = tt/R. 



S d GA k + (fc 2 + k 2 )B k = 0. 



(107) 



There will be non-trivial solutions only if the determi- 
nant of this system of equations is zero. This yields the 
dispersion relation 



SdGp 



k R T 



(108) 
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relating <Jk to the wavenumber k. The amplitudes Ak 
and Bk are determined by the initial condition. We see 
that Cfc is real so that the perturbation either grows or 
decays exponentially rapidly. The homogeneous phase 
will be spectrally stable if Uk < for all fc and it will be 
spectrally unstable if there exists one or several modes for 
which au > 0. We note that the dispersion relation (|108[) 
is the same as in an infinite domain [57J. However, in a 
finite domain, the allowed wavenumbers k are quantized 
while they are continuous in an infinite domain. 

According to Eq. (|108[) , the system will be unstable if 
there exists k ^ such that 



S d Gp k B T 



fc 2 + fc2 



(109) 



Therefore, a necessary condition of instability is that 

k B T* 



k B T S d Gp 



(110) 



where kf is the smallest non-zero wavenumber. For 
T > T* , the homogeneous distribution is stable for per- 
turbations with arbitrary wavenumbers. For T < T*, 
the homogeneous distribution is unstable for perturba- 
tions with wavenumbers 



2 S d Gmp 2 2 

For kg = 0, the critical temperature is 
k B T* S d Gp 



in 



kj 



and we recover the Jeans criterion 



k 2 < S d Gmp ^ fc2; 



k B T 



(111) 



(112) 



(113) 



In the general case, the instability criterion can be writ- 
ten 



k < kj k Q — k^ 



(114) 



We see that, for the screened Newtonian potential (ko ^ 
0), the instability occurs for larger wavelengths as com- 
pared to the Newtonian model (ko = 0). Let us introduce 
the notation 



kuT r = 



SdGmp 



k 2 



(115) 



which corresponds to the critical temperature in an infi- 
nite domain (kf =0). Since the dispersion relation (|108[) 
does not explicitly depend on kf, it is convenient to in- 
troduce the notation (|115[) . We have 



When T < T*, the system is unstable for the modes such 
that 



1 T~ l 



1/2 



k m (T). 



The growth rate can be written 

k B TP(k m (Tf-k 2 ) 



m 



fcg ~l~ k~ 



where 



1 T' 1 



1/2 



(117) 



(118) 



(119) 



It achieves its maximum value for k = k* (T) where 



fc*(T) = fc 



1/2 



1/2 



(120) 



The corresponding value of the growth rate is 



/m \ k B T c 2 
cr*(J ) = fe n 



■in 



1/2 



-I 2 



(121) 



The number of clusters that is expected to form in the 
linear regime is N(T) = Rj ' (2ir/k*(T)). For a fixed value 
of fen, this number increases as the temperature decreases. 
The behaviour of the different quantities defined above 
is represented in Figs. H9"l and ISTil 
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FIG. 49: Growth (a > 0) or decay (a < 0) rate as a function 
of the wavenumber k. The system is unstable for k < k m (T) 
and the maximum growth rate is reached for k — k*(T). The 
parameters have been scaled such that ko = 1, T c = 1, T = 
1/2. 



Let us consider some particular cases: 

• For T = T c , we have k m = 0, fc* = 0, c* = and 



l + (k f /k ) 



(116) 



Ok 



< 0. 



(122) 
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FIG. 50: Evolution of fc m , fc* and <r, as a function of the 
temperature. The parameters have been scaled such that ko = 
1 and T c = 1. 
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FIG. 51: Growth (a > 0) rate as a function of the wavenumber 
k in two limits: (i) the Newtonian limit fco = (and for 
which the maximum growth rate corresponds to k, — kf « 1 
(large scales), and (ii) the cold limit T = (and ko 7^ 0) for 
which the maximum growth rate corresponds to fe» — > +oo 
(small scales). 



Therefore, the system is stable. More generally, for T > 
T c , the system is stable. For T — > +oo. we have a k = 

m 

• For T = 0, we have k m — > +oo, fc» — > +oo, cr» — > 
kBTck^/m and 



k 2 k 2 



The growth rate is maximum for fc* — > +oo, i.e. for very 
small wavelengths A* — > 0. In that case, we expect a very 
large number of clusters in the linear regime. 

• For fcg = (modified Newtonian model) , we have 



a k = SdGp 



kuT, 



(124) 



The system is unstable for T < T* where the critical 
temperature is given by Eq. (|112l) . Furthermore, the 



unstable wavenumbers correspond to k < kj where the 
Jeans wavenumber is given by Eq. (|113[) . The growth 
rate is maximum for fc* = kf i.e. for the maximum 
wavelength A/ = 2ir/kf. In that case, we have only one 
cluster. The corresponding value of the growth rate is 
cr* = SdGp — kBTk'j/m. 

The two limit cases discussed above are illustrated in 

Fig. ED 



B. Thermodynamical stability 

We now analyze the thermodynamical stability of the 
homogeneous phase by using variational principles. Ba- 
sically, we have to solve the maximization problem J3]) 
in MCE and the minimization problem (|17p in CE. How- 
ever, for spatially homogeneous systems, it is shown in 
Appendix [X] that they are both equivalent to the mini- 
mization problem (f2"4"]l . Therefore, the system is stable iff 
the second order variations of free energy (|27|) are posi- 
tive definite for any perturbations Sp that conserve mass, 
i.e. J Spdr = 0. We are led therefore to considering the 



eigenvalue problem 



kgT 

5$x + Sp x = \5p x , 

pm 



ASQx - k 2 Q 6<5> x = S d G5 P x 



(125) 



(126) 



If all the eigenvalues A are positive, then the system is 
stable since S 2 F = | Xa\ > where the perturbation 
has been decomposed under the form Sp = J2\ a \3p\- If 
at least one eigenvalue is negative, the system is unstable 
since 5 2 F = ^A J {Sp\) 2 dr < for that perturbation. It 
is easy to see that the eigenfunctions are 

a r 1 a 

6p(r) = AMr), S*(r) = ^fc(r), (127) 

and that the corresponding eigenvalues are 



\(k) = - 



SdG 



fc 



pm ' 



(128) 



for all quantized k (see Sec. IV A[) . We note that 
fili k dr = -^ J AV'fe dr = -^§ ■ dS = 0, so that 
J Spdr = as required. Regrouping all these results, we 
conclude that the system is stable iff 



S d G 



k 2 



ko 



k B T 
pm 



< 0, 



(129) 



for all (quantized) k. This returns the stability condi- 
tion obtained in Sec. IV Al Therefore, the system is 
stable iff T > T*. If T < T*, the homogeneous phase 
is an unstable saddle point of free energy at fixed mass. 
This method proves the thermodynamical stability of the 
homogeneous phase in the canonical and microcanonical 
ensembles. This implies the stability with respect to the 
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mean field Kramers equation (f21j) , with respect to the 
Smoluchowski equation ([25)1 , with respect to the Keller- 
Segel model (|31|) and with respect to the kinetic equation 

(USD. 

We can now determine the values of the normalized 
inverse temperature 77* and normalized energy A* above 
which the homogeneous phase becomes unstable. Using 
Eqs. (jnj) and (|TTU|) . we get 



1 



(kj 



kl)R\ 



(130) 



Wc obtain 



7T 2 + M 2 = 9.8696044 +/i 2 (d = 1), (131) 



V* c =\{jli+V 2 ) = 7.3410008+^ 



Vc = \^\+^) = 6.7302445+ ^ 



(d = 2),(132) 



(d = 3). (133) 



The corresponding critical energy is given by Eq. (|63[) 
for the modified Newtonian model and by ([57]) for the 
screened Newtonian model. 



VI. CONCLUSION 

In this paper, we have completed the description of 
phase transitions in self-gravitating systems and bacte- 
rial populations. We have introduced generalized models 
in which the ordinary Poisson equation is modified so 
as to allow for the existence of a spatially homogeneous 
phase. This avoids the Jeans swindle and leads to a great 
variety of microcanonical and canonical phase transitions 
between homogeneous and inhomogeneous states. These 
generalized models can have application in chemotaxis 
where the degradation of the chemical leads to a shield- 
ing of the interaction and in cosmology where the ex- 
pansion of the universe creates a sort of "neutralizing 
background" . In this paper, we have only considered 
equilibrium states. In future works, we shall study the 
dynamics of some simple models for which the present 
study can be a useful guide. 

Our study also allows to explore the link between cos- 
mology where one studies the evolution of the universe 
as a whole [l2j and stellar dynamics where one studies 
the structure of individual galaxies [59| . The description 
of phase transitions in these two disciplines is usually 
very different [60j . However, our study allows to make 
some basic connections. In cosmology, one usually starts 
from an infinite homogeneous distribution and study the 
appearance of clusters representing galaxies. Our ther- 
modynamical approach shows that, indeed, the homo- 
geneous phase is unstable for sufficiently low tempera- 
tures and energies and leads to clusters. The formation 
of these clusters can be studied by making a linear sta- 
bility analysis of the Vlasov or Euler equations. Then, in 



the nonlinear regime, the system is expected to achieve 
a statistical equilibrium state due to violent relaxation 
or collisional relaxation (finite N effects) This cor- 
responds to the inhomogeneous phase. In d = 1, there 
exists an equilibrium state for any value of energy and 
temperature. For low energies and temperatures, it is 
spatially inhomogeneous. In the core of the cluster, the 
density is so high that we can disregard the effect of 
the neutralizing background. In that case, the statistical 
equilibrium state (representing a "galaxy") is described 
by the Camm solution like in ID stellar dynamics. In 
d = 3, there is no inhomogeneous equilibrium state and, 
for sufficiently small energies and temperatures, the sys- 
tem undergoes a gravothermal catastrophe or an isother- 
mal collapse. In d = 2, the situation is intermediate. 
There exists an equilibrium state in the microcanonical 
ensemble for all energies while in the canonical ensemble 
no equilibrium state exists at low temperatures. Similar 
behaviours occur in chemotaxis and will be investigated 
in future papers. Note that for self-gravitating systems, 
the proper statistical ensemble is the microcanonical en- 
semble while in chemotaxis (or for the academic model of 
self-gravitating Brownian particles) the proper statistical 
ensemble is the canonical one. It is therefore interesting 
to study these two systems in parallel to describe the 
analogies and differences between statistical ensembles. 



Appendix A: Equivalent but simpler optimization 
problems 

In this Appendix, following the approach of Padman- 
abhan [22j and Chavanis [27]], we shall reduce the op- 
timization problems ([3]) and (|17[) to simpler forms. In 
particular, we shall show that these optimization prob- 
lems for /(r, v) are equivalent to optimization problems 
for p(r). 



1. Microcanonical ensemble 

To solve the maximization problem ([3]) we can pro- 
ceed in two steps. Wc first maximize the entropy at fixed 
energy, mass and density profile p(r). Since the spec- 
ification of p(r) determines the mass and the potential 
energy, this is equivalent to maximizing the entropy at 
fixed kinetic energy and density profile. Writing 

5S-/35 (^j f^-drdv^j - J X(r)s(J fdv^j dr = 0, 



(Al) 



this leads to the Maxwellian distribution function 



/(r,v) 



/ m 



d/2 



p(r)e a B r , (A2) 



\2nk B T , 

which is the global entropy maximum with the previ- 
ous constraints since 5 2 S = — J ^fln ^ rc ^ v < (the con- 
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straints are linear in / so that their second variations 
vanish). We can now express the mass, the entropy and 
the energy in terms of p(r) and T. Up to unimportant 
constants, we obtain 



5" = -Nk B InT 



^In^dr, 



(A3) 



E=^Nk B T + ± J p(v,t)u(r,r')p(v',t)drdr' 

pV dr. (A4) 



We now have to solve the maximization problem 

max {S[p\ | E[p] = E, M[p] = AI} . (A5) 

p 

Finally, the solution of ([3]) is given by the distribution 
function (|A2[) with the density profile that is solution 
of l|A5p . Let us compute the variations of entropy and 
energy up to second order. We have 

AS= W£ -k B [(lnt + i yJL dr 
2 1 J \ to / to 

"-k B [^dr, (A6) 
J 2pm 



d ,r, (ST 
-Nk B — 
4 V T 



AE = ^Nk B 5T + J $5pdr+^J 5p5<Pdr. (A7) 

Using the conservation of energy AE = to eliminate 
5T, we obtain 

AS = [ $5pdr -k B [ (hi + l) ^ dr 
T J J \ m / to 

-^I spmdr -dd^(I Mpdr ) 

-k B [^dr. (A8) 
J 2pm 

Let us consider the first order variations. At first order, 
we have 

SS = ~ [ <$>5pdr -k B [ fin + l) ^ dr. (A9) 

T J J V TO / TO 

The conservation of mass can be taken into account by in- 
troducing a Lagrange multiplier. Writing the variational 
principle as 

5S - aSM = 0, (A10) 
we obtain the mean field Boltzmann distribution 

p = A'e~^, (All) 



where $(r) is given by Eq. ([7]). Combining Eq. (|A11|) 
with Eq. (| A2|) . we recover the mean field Maxwell- 
Boltzmann distribution (|10[) . However, the present ap- 
proach allows us to simplify the condition of thcrmody- 
namical stability. Indeed, the system is stable in the 



microcanonical ensemble iff the second order variations 
of entropy (|A8p are negative definite 



(Spf 
2 pm 
1 



1 f 

dr — — / 8p5Q dr 



dNk B T 2 



<$>Spdr) <0 



(A12) 



for any perturbation Sp that conserves mass at first order, 
i.e. J Spdr = (the conservation of energy has automati- 
cally been taken into account in the previous derivation) . 
This stability criterion is equivalent to the stability crite- 
rion (|12|) but it is simpler because it is expressed in terms 
of the density instead of the distribution function. 



2. Canonical ensemble 

To solve the maximization problem (|17p we can pro- 
ceed in two steps. We first minimize the free energy at 
fixed mass and density profile p(r). This is equivalent to 
minimizing the free energy at fixed density profile. Writ- 
ing 

SF + T J \{r)6 (^J fdv^j dr = 0, (A13) 



this leads to the Maxwcllian distribution function 

d/2 2 



/(r,v) 



/ TO 



V 2-nk B T 



p(r) e 2k 



(A14) 



which is the global minimum of free energy with the pre- 
vious constraint since S 2 F = T J ^Jm ^ rc ' v > (^ ne 
constraints are linear in / so that their second variations 
vanish). We can now express the free energy in terms of 
p(r). Up to unimportant constants, we get 

F=\J p(v,t)u(r,r')p(v',t) drdr' 

pV dr + k B T ( — In — dr. (A15) 
J m m 

We now have to solve the minimization problem 

min {F[p\ | M[p] = M} . (A16) 
p 

Finally, the solution of (|17p is given by the distribution 
function (|A14|) with the density profile that is solution of 
(|A16p . The first variations 

5F + aTSAI = 0, (A17) 

lead to the mean field Boltzmann distribution 

p = A'e~^, (A18) 



where $(r) is given by Eq. ([7]). Combining Eq. (|A18[) 
with Eq. (|AT4|) , we recover the mean field Maxwell- 
Boltzmann distribution (|10p . However, the present ap- 
proach allows us to simplify the condition of thcrmody- 
namical stability. Indeed, the system is stable in the 
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canonical ensemble iff the second order variations of free 
energy (|A15p are positive definite 



2 I 



k B T f(5p) 



2p 



dr > 0, 



(A19) 



for any perturbation dp that conserves mass at first or- 
der, i.e. J Spdr = 0. This stability criterion is equivalent 
to the stability criterion (f2"U| but it is simpler because it 
is expressed in terms of the density instead of the distri- 
bution function. 

Remark: From the stability criteria (|A12[) and (|A19|) , 
we clearly see that canonical stability implies micro- 
canonical stability (but not the converse). Indeed, since 
the last term in Eq. (|A12[) is negative, it is clear that 
if inequality (|A19|) is satisfied then inequality (|A12[) is 
automatically satisfied. In general, this is not reciprocal 
and we may have ensembles inequivalence. However, if 
we consider a spatially homogeneous system for which $ 
is uniform, the last term in Eq. (|A12[) vanishes (since 
the mass is conserved) and the stability criteria (|A12[) 
and (|A19[) coincide. Therefore, for spatially homoge- 
neous systems, we have ensembles equivalence. 



Appendix B: Explicit expressions of the potential 

In this Appendix, we limit ourselves to the case d — 1 
although the results can be easily generalized to any 
dimension. Using standard methods, we can obtain 
the Green function associated with the screened Poisson 
equation (|64[) in a box with Neumann boundary condi- 
tions. Then, we find that the potential is explicitly given 
by 



$(x) = / p(x')u{x,x')dx' 

-R 



(Bl) 



with 



u(x, x 1 ) 



G 



ko sinh(2fcoi?) 
x (cosh [k {2R -\x- x'\)] + cosh [k (x + x')}) 

In an infinite domain (R — > +oo), we obtain 



(B2) 



G 

u(\x — x'\) = — — e 

Kq 



— ko\x — x'\ 



(B3) 



Similarly, for the modified Newtonian model (|4"3"]l with 
Neumann boundary conditions, the potential is explicitly 
given by 



${x) = G (p-p){x')\x-x'\dx', (B4) 

J-R 

and this expression remains valid in an infinite domain 
(R +oo). 



Appendix C: External potential for the modified 
Newtonian model 



For the modified Newtonian model (|43j) . the potential 
can be written 

*(r) = J (p-p)(r')u(r,v')dr', (CI) 

where u(r,r') is the Green function of the Poisson equa- 
tion with Neumann boundary conditions. Comparing Eq. 
(|C1[) with Eq. ([7]), we find that the external potential is 



V(r) = -Pj u(r,r')dr'. (C2) 

Using Eq. (|C2[) . the potential energy can be written 

W =\j P^dr-^p J p(r)u(r, r') drdr' . (C3) 

Interchanging the dummy variables r and r' and using 
the symmetry u(r',r) = u(r,r'), we get 

w = \ J dr ~ \p J P( r ') u (r, r ') drdr'. (C4) 

Finally, using Eq. (O, we obtain 

W= \J(p-p)®dr+^pJv(r)dr. (C5) 

Therefore, the potential energy is given by Eq. (|53|) 
up to an unimportant additive constant \fp J V(r) dr = 
— \~p 2 j u(r,r')drdr' . 

In d = 1, according to Eq. (|B4|) . the potential of inter- 
action is u = G\x — x'\. Therefore, the external potential 
is explicitly given by 

V(x) = ~pG{x 2 +R 2 ). (C6) 

The additive constant in the energy is 



l --p J ^ V(x)dx = -^Gp 2 R s . 



(C7) 



Appendix D: The minimum energy 

Let us consider the modified Newtonian model (|43|) in 
d = 1. At T = 0, the density profile is a Dirac peak 
p = AIS(x) and the energy ([56]) is 



1 f B fd$\ 2 , 



(Dl) 



For a symmetric density profile, the modified Poisson 
equation can be integrated into 

$'(x) = 2G / p{x) dx - 2Gpx, (D2) 
Jo 
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which is the appropriate Gauss theorem. If all the mass 
is concentrated at x = 0, we obtain 



$'(x) = GM (sign(x) 



R 



(D3) 



Substituting this expression in Eq. (|D1[) . we obtain E 
-GM 2 R/6. The total normalized energy is therefore 



A« =-. 

max g 



(D4) 



This corresponds to the minimum energy of the branch 
n = 1. 

Let us consider the screened Newtonian model (|64|) 
in d = 1. At T = 0, the density profile is a Dirac peak 
p = AIS(x) and the energy (J7¥J) w E = ^M$ - According 
to Eq. (|B2[) , the potential created in i by a mass M 
located at x' = is 



$(x) = - 



GM cosh [fc (-R- M 



fc 



sinh(fc i?) 



(D5) 



where we have used elementary trigonometric identi- 
ties to simplify the expression. This leads to $o = 
-GM/(fc tanh(fc i?)) and E = -GM 2 /{2k tanh(fco-R)). 
The total normalized energy is therefore 



A« = - . 

max 2/itanh(p) 



(D6) 



This corresponds to the minimum energy of the branch 
n = 1. 



Appendix E: Approximate expressions of the density 
profile 

In d — 1, the screened Emden equation can be 
written 



with 



d 2 i> 



- A + K 2 ip = - 



dV 

dip' 



2„/,2 



(El) 



(E2) 



This is similar to the equation of motion of a particle 
of unit mass in a potential V^) where ip plays the role 
of position and £ the role of time. Using the boundary 
condition ip = ip' = at £ = 0, we find that the first 
integral (pseudo energy) is 



(E3) 



This first order differential equation can be easily inte- 
grated until £ = a, which formally solves the problem. 

Let us consider the limit po ~ ^ +°° corresponding to 
a — > +oo. In the inner region, the term dominates 



and Eq. (|E1[) reduces to the ordinary Emden equation 
whose solution is the Camm profile [29|, [6l[ : 



-i> - 



1 



cosh 2 (f/V2)' 



(E4) 



In the outer region, the term e ^ can be neglected and 
Eqs. ((El) and (jE3j) reduce to 



d 2 ip 

He 



X + K 2 lp, 



1 [dip 



2\d4 



1 



2„/.2 



+ A?/> — — K Ip 



1. 



(E5) 



(E6) 



The boundary condition at the wall is ip'(a) = 0. Substi- 
tuting this result in Eq. (|E6|) . we get \ip(a) — \n 2 ip{a) 2 = 
1. The physical solution of this equation is ^(a) 
: A - \/A 2 - 2k 2 ) /k 2 . Solving Eq. (|E5|) with these bound- 
ary conditions, we find that 



A 



\ yj\ 2 - 2k 2 cosh [k{£ - a)] . (E7) 



, 2 



The matching of the outer solution with the inner solu- 
tion implies that ip outer (fi) = 0. Using na = /i, we obtain 



A V2~n 1 ( a -r+oo) 
tanh(p) a ' 



(E8) 



Finally, substituting the inner profile (jE4J) in Eq. (|73[) . 
we obtain at leading order 



^2a, (a — > +oo). 



(E9) 



For a — > +oo, the density profile tends to a Dirac 
peak p = M8(x). The potential energy reduces to 
W = ±M$ . Using Eqs. (JTHJ), CE3) and « = /x/a, we 
recover Eq. (|D6[) . 

The modified Emden equation (|47|) can be studied sim- 
ilarly. In fact, most of the preceding results remain valid 



by taking k = 0. The potential is V(ip) 



Xip. 



It has a minimum at tpo = — In A so that the solution 
of the Emden equation (with energy £ = 1) oscil- 
lates around this value. Integrating Eq. (|E3|) . the density 
profiles of the solutions of the branch n = 1 are given by 



<7.r 



(E10) 



with £, < a. The half-period of the oscillations of the 
function "0(C) is 



L 



i/i(a) 



dx 



^/2(1 - e"* - Ax) ' 



(Ell) 



where VK ) is solution of e + Xip(a) = 1 obtained 
from Eq. (|E3I) with ip'{a) = 0. Let us now consider 
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the limit p$ — > +00. The inner solution is given by the 
Camm profile (|E4[) and the outer solution is 

m = \-\Kt-<*) 2 , (E12) 

which is consistent with Eq. (IE7I) when k — > 0. The 
matching condition ipouter (0) = then yields 




a 



Using Eq. ([52]) . we obtain at leading order 

r\ ~ V2a, {a +00). (E14) 



with the boundary conditions • n = on the bound- 
ary. The boundary conditions determine the allowable 
wavenumbers k 2 = S d G(3mp. They take discrete val- 
ues k = k„ (see Sec. |V| which in turn determine dis- 
crete values of the temperature T n . The first point of 
bifurcation corresponds to the smallest wavenumber kf. 
This is associated with the critical temperature (|112p at 
which the homogeneous branch becomes unstable. Other 
branches of bifurcations appear at smaller temperatures. 
They correspond to successive quantized values k n of the 
wavenumber. 

For the screened Newtonian model, the differential 
equation determining the field $(r) at statistical equi- 
librium can be written 



Appendix F: The bifurcation point 

In this Appendix, we shall determine the point at 
which the spatially homogeneous branch bifurcates to the 
spatially inhomogeneous branch and show that it coin- 
cides with the point at which the spatially homogeneous 
branch becomes unstable (see Sec. |V]). For a more de- 
tailed theory of bifurcations, we refer to the paper of 
Schaff H. 

For the modified Newtonian model, the differential 
equation determining the field $(r) at statistical equi- 
librium can be written 

A$ = S d G (Ae' " 1 * - p) . (Fl) 

The homogeneous solution corresponds to p = ~p, $ = 
and A — p. Considering a small perturbation $ = + 
cj)(r) with (/) <C 1 around the homogeneous solution and 
linearizing the differential equation (|F1[) . we obtain 

A0 + S d G/3mp4> = 0, (F2) 



A$ - k 2 ^ = S d GAe~' )m *. (F3) 



The homogeneous solution corresponds to p = const., 
<3> = const, with — /jq$ = S d Gp. Considering a small per- 
turbation $ = const. +<f>(r) with (j> <C 1 around the homo- 
geneous solution and linearizing the differential equation 
(|F3j) . we obtain 

A(j) + {S d Gf3mp - kl)(j) = 0, (F4) 
with the boundary conditions • n = on the bound- 
ary. The boundary conditions determine the allowable 
wavenumbers k 2 = S d G(3mp — k 2 . The first point of bi- 
furcation corresponds to the smallest wavenumber kf (see 
Sec. |V]) . This is associated with the critical temperature 
(|110|) at which the homogeneous branch becomes unsta- 
ble. Other branches of bifurcations appear at smaller 
temperatures. They correspond to the successive quan- 
tized values k n of the wavenumber. 
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